From e6abb2226d1681a4dd661a0b2f018392c1e7ffb6 Mon Sep 17 00:00:00 2001 From: zh-qs Date: Thu, 19 Jan 2023 00:47:51 +0100 Subject: [PATCH 01/13] add performance tests --- engine/test/AdvancedIntegral.cu | 182 +++++++++++++++++++++++++++----- 1 file changed, 155 insertions(+), 27 deletions(-) diff --git a/engine/test/AdvancedIntegral.cu b/engine/test/AdvancedIntegral.cu index 7968fb0..0a4cc96 100644 --- a/engine/test/AdvancedIntegral.cu +++ b/engine/test/AdvancedIntegral.cu @@ -1,9 +1,114 @@ #include "IntegralCommons.cuh" -#define ADVANCED_INTEGRAL_TEST(_name, _integral, _expected_result) \ +#include +#include +#include +#include +#include + +#include "Evaluation/Integrator.cuh" +#include "Parser/Parser.cuh" +#include "Simplify.cuh" +#include "Symbol/Macros.cuh" +#include "Symbol/Symbol.cuh" + +#define ADVANCED_INTEGRAL_TEST_CUDA_ONLY(_name, _integral, _expected_result) \ INTEGRAL_TEST(AdvancedIntegral, _name, _integral, _expected_result) +#define ADVANCED_INTEGRAL_TEST_TIME(_name, _integral, _expected_result) \ + TEST_F(AdvancedIntegral, _name) { test_integral_with_external_tools(_integral); } + +#define ADVANCED_INTEGRAL_TEST(_name, _integral, _expected_result) \ + ADVANCED_INTEGRAL_TEST_TIME(_name, _integral, _expected_result) + namespace Test { + namespace { + std::string exec_and_read_output(const std::string cmd) { + // sciagniete z SolverProcessManager.cu + std::array buffer{}; + std::string result; + std::unique_ptr pipe(popen(cmd.c_str(), "r"), pclose); + if (!pipe) { + throw std::runtime_error("popen() failed!"); + } + while (fgets(buffer.data(), buffer.size(), pipe.get()) != nullptr) { + result += buffer.data(); + } + return result; + } + + void capitalize_substring(std::string& string, const std::string& what) { + if (string.size() < what.size()) { + return; + } + for (size_t i = 0; i <= string.size() - what.size(); ++i) { + if (string.substr(i, what.size()) == what && + (i == 0 || !std::isalpha(string[i - 1])) && + (i == string.size() - what.size() || !std::isalpha(string[i + what.size()]))) { + string[i] = string[i] - 'a' + 'A'; + } + } + } + + void replace(std::string& string, const std::string& from, const std::string& to) { + if (string.size() < from.size()) { + return; + } + for (ssize_t i = 0; i <= string.size() - from.size(); ++i) { + if (string.substr(i, from.size()) == from) { + string = string.substr(0, i) + to + string.substr(i + from.size()); + } + } + } + + std::string make_mathematica_command(const std::string& integral) { + // /media/cmonbrug/DATA/mathematica/Executables/wolframscript -code + // 'expr=ToExpression["sin(Pi)+x",TraditionalForm];{t,b}=AbsoluteTiming[Integrate[expr,x]];t' + std::string expression = integral; + capitalize_substring(expression, "e"); + capitalize_substring(expression, "pi"); + capitalize_substring(expression, "sqrt"); + return fmt::format( + R"(wolframscript -code 'expr=ToExpression["{}",TraditionalForm];{{t,b}}=AbsoluteTiming[Integrate[expr,x]];t')", + expression); + } + + std::string make_sympy_command(const std::string& integral) { + std::string expression = integral; + capitalize_substring(expression, "e"); + replace(expression, "^", "**"); + replace(expression, "arctan", "atan"); + replace(expression, "arccos", "acos"); + replace(expression, "arcsin", "asin"); + replace(expression, "sgn", "sign"); + //fmt::print(R"(Call: python3 -c 'from sympy import *;x=Symbol("x");print(utilities.timeutils.timed(lambda:integrate({},x))[1])')" "\n", + // expression); + return fmt::format( + R"(python3 -c 'from sympy import *;x=Symbol("x");t=Symbol("t");print(utilities.timeutils.timed(lambda:integrate({},x))[1])')", + expression); + } + + void test_integral_with_external_tools(const std::string& integral_str) { + const auto integral = Sym::integral(Parser::parse_function(integral_str)); + + Sym::Integrator integrator; + + const clock_t start = clock(); + const auto result = integrator.solve_integral(integral); + const clock_t end = clock(); + + const double cusymint_seconds = static_cast(end - start) / CLOCKS_PER_SEC; + + auto mathematica_result = exec_and_read_output(make_mathematica_command(integral_str)); + auto sympy_result = exec_and_read_output(make_sympy_command(integral[1].to_string())); + + printf("Cusymint time: %f\n" + "Mathematica time: %s\n" + "Sympy time: %s\n", + cusymint_seconds, mathematica_result.c_str(), sympy_result.c_str()); + } + } + class AdvancedIntegral : public IntegrationFixture {}; // easy excercises from http://www.math.uni.wroc.pl/~ikrol/zbiorek7.pdf @@ -19,24 +124,29 @@ namespace Test { ADVANCED_INTEGRAL_TEST(CosineByParts, "x cos(x)", "x sin(x)+cos(x)") ADVANCED_INTEGRAL_TEST(ExponentByParts, "x e^x", "x e^x-e^x") ADVANCED_INTEGRAL_TEST(ExponentByParts2, "x^2e^x", "e^x(x^2-2x+2)") - //ADVANCED_INTEGRAL_TEST(NegativeExponentByParts, "x e^(-x)", "-e^(-x)(x+1)") // linear substitution for nested-vars-only required + // ADVANCED_INTEGRAL_TEST(NegativeExponentByParts, "x e^(-x)", "-e^(-x)(x+1)") // linear + // substitution for nested-vars-only required ADVANCED_INTEGRAL_TEST(CosineByParts2, "x^2cos(x)", "x^2sin(x)+2x cos(x)-2sin(x)") ADVANCED_INTEGRAL_TEST(LogarithmTimesRoot, "sqrt(x)ln(x)", "2/3x^(3/2)(ln(x)-2/3)") ADVANCED_INTEGRAL_TEST(LogarithmDividedByPower, "ln(x)/x^4", "-1/(3x^3)(ln(x)+1/3)") ADVANCED_INTEGRAL_TEST(VeryLongPolynomial, "(2x-1)^20", "1/42(2x-1)^21+1/42") ADVANCED_INTEGRAL_TEST(CosineOfLinearFunction, "cos(3x-1)", "sin(3x-1)/3") - //ADVANCED_INTEGRAL_TEST(RationalWithSubstitution, "2x/(5x^2+1)", "1/5ln(5x^2+1)") // substitution t=x^2+a required + // ADVANCED_INTEGRAL_TEST(RationalWithSubstitution, "2x/(5x^2+1)", "1/5ln(5x^2+1)") // + // substitution t=x^2+a required ADVANCED_INTEGRAL_TEST(SimpleRationalWithExponent, "e^x/(3+e^x)", "ln(3+e^x)") - //ADVANCED_INTEGRAL_TEST(ExponentWithSquare, "x*e^(-x^2)", "-1/2e^(-x^2)") // substitution t=x^2+a required - //ADVANCED_INTEGRAL_TEST(ExponentOfInverse, "e^(1/x)/x^2", "-e^(1/x)") // substitution t=ax^b+c required (dx = 1/b((t-c)/a)^((1-b)/b)) + // ADVANCED_INTEGRAL_TEST(ExponentWithSquare, "x*e^(-x^2)", "-1/2e^(-x^2)") // substitution + // t=x^2+a required ADVANCED_INTEGRAL_TEST(ExponentOfInverse, "e^(1/x)/x^2", "-e^(1/x)") // + // substitution t=ax^b+c required (dx = 1/b((t-c)/a)^((1-b)/b)) ADVANCED_INTEGRAL_TEST(Tangent, "tg(x)", "-ln(cos(x))") ADVANCED_INTEGRAL_TEST(TangentDividedByCosSquared, "tg(x)/cos^2(x)", "tg^2(x)/2") ADVANCED_INTEGRAL_TEST(CosTimesExpSine, "cos(x)e^sin(x)", "e^sin(x)") ADVANCED_INTEGRAL_TEST(CosineDividedByRootOfLinSine, "cos(x)/sqrt(1+sin(x))", "2sqrt(1+sin(x))") - //ADVANCED_INTEGRAL_TEST(TangentDerivativeWithSubstitution, "x^3/cos^2(x^4)", "tg(x^4)/4") // substitution t=ax^b+c required (dx = 1/b((t-c)/a)^((1-b)/b)) + // ADVANCED_INTEGRAL_TEST(TangentDerivativeWithSubstitution, "x^3/cos^2(x^4)", "tg(x^4)/4") // + // substitution t=ax^b+c required (dx = 1/b((t-c)/a)^((1-b)/b)) ADVANCED_INTEGRAL_TEST(ExponentOfLinear, "6^(1-x)", "-6^(1-x)/ln(6)") - //ADVANCED_INTEGRAL_TEST(ArctanSubtitution, "1/(1+x^2)/arctan(x)", "ln(arctan(x))") // substitution t=arctan(x) required + // ADVANCED_INTEGRAL_TEST(ArctanSubtitution, "1/(1+x^2)/arctan(x)", "ln(arctan(x))") // + // substitution t=arctan(x) required // task 4 contains examples which require nested integrals (not implemeneted) @@ -45,39 +155,57 @@ namespace Test { ADVANCED_INTEGRAL_TEST(SimpleRationalWithSquare, "(sqrt(x)-2)^2/x^2", "-4/x+8/sqrt(x)+ln(x)") ADVANCED_INTEGRAL_TEST(SineCosineSquaredInDenominator, "1/(sin^2(x)cos^2(x))", "-1/tg(x)+tg(x)") - //ADVANCED_INTEGRAL_TEST(LogarithmDividedByX, "ln^5(x)/x", "ln^6(x)/6") // substitution t=ln(x) or cyclic integrals required - //ADVANCED_INTEGRAL_TEST(ArctangentWithX3Substitution, "x^3/(x^8+1)", "arctg(x^4)/4") // substitution t=x^4 required - //ADVANCED_INTEGRAL_TEST(LogarithmMultipliedByX, "x^n*ln(x)", "x^(n+1)ln(x)/(n+1)-x^(n+1)/(n+1)^2") // non-numeric powers integration by parts required + // ADVANCED_INTEGRAL_TEST(LogarithmDividedByX, "ln^5(x)/x", "ln^6(x)/6") // substitution t=ln(x) + // or cyclic integrals required ADVANCED_INTEGRAL_TEST(ArctangentWithX3Substitution, + // "x^3/(x^8+1)", "arctg(x^4)/4") // substitution t=x^4 required + // ADVANCED_INTEGRAL_TEST(LogarithmMultipliedByX, "x^n*ln(x)", + // "x^(n+1)ln(x)/(n+1)-x^(n+1)/(n+1)^2") // non-numeric powers integration by parts required ADVANCED_INTEGRAL_TEST(ArctangentMultipliedByX, "x*arctg(x)", "x^2/2*arctg(x)-x/2+arctg(x)/2") - //ADVANCED_INTEGRAL_TEST(Arcsine, "arcsin(x)", "x*arcsin(x)+2*sqrt(1-x^2)") // substitution t=x^2+a required - //ADVANCED_INTEGRAL_TEST(LogarithmWithX2Substitution, "x*ln(x^2+1)", "(x^2+1)ln(x^2+1)-x^2-1") // substitution t=x^2+a required - //ADVANCED_INTEGRAL_TEST(CyclicIntegral1, "e^x*sin(2x)", "e^x*sin(2x)/5 - e^x*2cos(2x)/5") // cyclic integrals required - //ADVANCED_INTEGRAL_TEST(CyclicIntegral2, "sin(ln(x))", "-x/2(cos(ln(x))-sin(ln(x)))") // cyclic integrals required - //ADVANCED_INTEGRAL_TEST(CotangentDividedByLogOfSine, "ctg(x)/ln(sin(x))", "ln(ln(sin(x)))") // substitution t=ln(x) required - //ADVANCED_INTEGRAL_TEST(InvertedCosH, "1/(e^x+e^-x)", "arctg(e^x)") // extract_function required + // ADVANCED_INTEGRAL_TEST(Arcsine, "arcsin(x)", "x*arcsin(x)+2*sqrt(1-x^2)") // substitution + // t=x^2+a required ADVANCED_INTEGRAL_TEST(LogarithmWithX2Substitution, "x*ln(x^2+1)", + // "(x^2+1)ln(x^2+1)-x^2-1") // substitution t=x^2+a required + // ADVANCED_INTEGRAL_TEST(CyclicIntegral1, "e^x*sin(2x)", "e^x*sin(2x)/5 - e^x*2cos(2x)/5") // + // cyclic integrals required ADVANCED_INTEGRAL_TEST(CyclicIntegral2, "sin(ln(x))", + // "-x/2(cos(ln(x))-sin(ln(x)))") // cyclic integrals required + // ADVANCED_INTEGRAL_TEST(CotangentDividedByLogOfSine, "ctg(x)/ln(sin(x))", "ln(ln(sin(x)))") // + // substitution t=ln(x) required ADVANCED_INTEGRAL_TEST(InvertedCosH, "1/(e^x+e^-x)", + // "arctg(e^x)") // extract_function required ADVANCED_INTEGRAL_TEST(Absolute, "abs(x)", "x*abs(x)/2") // abs required - ADVANCED_INTEGRAL_TEST(MaxOfOneAndSquare, "(abs(1-x^2)+1+x^2)/2", "(sgn(1-x^2)(x-x^3/3)+x+x^3/3)/2") // abs required + ADVANCED_INTEGRAL_TEST(MaxOfOneAndSquare, "(abs(1-x^2)+1+x^2)/2", + "(sgn(1-x^2)(x-x^3/3)+x+x^3/3)/2") // abs required // tasks 5-7 involve rational integrals // https://pages.mini.pw.edu.pl/~dembinskaa/www/?download=Inf_I_PowtKol3_2022-2023.pdf, task 11 ADVANCED_INTEGRAL_TEST(TangentSquared, "tg^2(x)", "tg(x)-arctg(tg(x))") // may sometimes fail ADVANCED_INTEGRAL_TEST(LongPolynomial, "(2x-3)^10", "1/22(2x-3)^11+3^11/22") - //ADVANCED_INTEGRAL_TEST(SquareRootInDenominator, "1/(2+sqrt(x))", "2sqrt(x)-4ln(sqrt(x)+2)") // substitution t=sqrt(x)+a required - //ADVANCED_INTEGRAL_TEST(InvertedCosH2, "1/cosh(x)", "2arctg(e^x)") // extract_funcion required + // ADVANCED_INTEGRAL_TEST(SquareRootInDenominator, "1/(2+sqrt(x))", "2sqrt(x)-4ln(sqrt(x)+2)") + // // substitution t=sqrt(x)+a required ADVANCED_INTEGRAL_TEST(InvertedCosH2, "1/cosh(x)", + // "2arctg(e^x)") // extract_funcion required ADVANCED_INTEGRAL_TEST(Sine5Cosine, "sin^5(x)cos(x)", "sin^6(x)/6") - //ADVANCED_INTEGRAL_TEST(ExpressionWithSquareInDenominator, "x/(x^2-1)^(3/2)", "-1/sqrt(x^2-1)") // substitution t=ax^b+c required - //ADVANCED_INTEGRAL_TEST(ArcsineWithDerivative, "arcsin^2(x)/sqrt(1-x^2)", "arcsin^3(x)/3") // substitution t=arcsin(x) required - //ADVANCED_INTEGRAL_TEST(RootInDenominator, "1/(x^(1/3)+1)", "3/2*x^(2/3)-3x^(1/3)+3ln(x^(1/3)+1)") // substitution t=ax^b+c required (dx = 1/b((t-c)/a)^((1-b)/b)) + // ADVANCED_INTEGRAL_TEST(ExpressionWithSquareInDenominator, "x/(x^2-1)^(3/2)", + // "-1/sqrt(x^2-1)") // substitution t=ax^b+c required + // ADVANCED_INTEGRAL_TEST(ArcsineWithDerivative, "arcsin^2(x)/sqrt(1-x^2)", "arcsin^3(x)/3") // + // substitution t=arcsin(x) required ADVANCED_INTEGRAL_TEST(RootInDenominator, "1/(x^(1/3)+1)", + // "3/2*x^(2/3)-3x^(1/3)+3ln(x^(1/3)+1)") // substitution t=ax^b+c required (dx = + // 1/b((t-c)/a)^((1-b)/b)) ADVANCED_INTEGRAL_TEST(Logarithm, "ln(x)", "x ln(x) - x") ADVANCED_INTEGRAL_TEST(XTimesCosine, "x*cos(x)", "x*sin(x)+cos(x)") - //ADVANCED_INTEGRAL_TEST(XTimesExponential, "x^2e^(1-x)", "-(x^2+2x+2)e^(1-x)") // linear substitution for nested-vars-only required - //ADVANCED_INTEGRAL_TEST(SquareTimesExponential, "8x^2e^(4-x^3)", "-8/3e^(4-x^3)") // substitution t=x^2+a required + // ADVANCED_INTEGRAL_TEST(XTimesExponential, "x^2e^(1-x)", "-(x^2+2x+2)e^(1-x)") // linear + // substitution for nested-vars-only required ADVANCED_INTEGRAL_TEST(SquareTimesExponential, + // "8x^2e^(4-x^3)", "-8/3e^(4-x^3)") // substitution t=x^2+a required - //ADVANCED_INTEGRAL_TEST(SineSquared, "sin^2(x)", "x/2-sin(2x)/4") // rational integrals required + // ADVANCED_INTEGRAL_TEST(SineSquared, "sin^2(x)", "x/2-sin(2x)/4") // rational integrals + // required ADVANCED_INTEGRAL_TEST(Sine5, "sin^5(x)", "-cos^5(x)/5+2/3cos^3(x)-cos(x)") ADVANCED_INTEGRAL_TEST(Sine4Cos3, "sin^4(x)cos^3(x)", "sin^5(x)/5-sin^7(x)/7") - //ADVANCED_INTEGRAL_TEST(Sine4Cos2, "sin^4(x)cos^2(x)", "1/6sin^5(x)cos(x)+1/24sin^3(x)cos(x)-1/32sin(x)cos(x)+1/16x") // rational integrals required - //ADVANCED_INTEGRAL_TEST(Cos3XCos5X, "cos(3x)cos(5x)", "1/16sin(8x)+1/4sin(2x)") // more trigonometric identitied required + // ADVANCED_INTEGRAL_TEST(Sine4Cos2, "sin^4(x)cos^2(x)", + // "1/6sin^5(x)cos(x)+1/24sin^3(x)cos(x)-1/32sin(x)cos(x)+1/16x") // rational integrals required + // ADVANCED_INTEGRAL_TEST(Cos3XCos5X, "cos(3x)cos(5x)", "1/16sin(8x)+1/4sin(2x)") // more + // trigonometric identitied required + + // own tests + ADVANCED_INTEGRAL_TEST(ETower, "e^x*e^e^x*e^e^e^x*e^e^e^e^x*e^e^e^e^e^x*e^e^e^e^e^e^x", + "e^e^e^e^e^e^x") } \ No newline at end of file From 34d2b9d2afce79f4e4caca0157be444e6d938f52 Mon Sep 17 00:00:00 2001 From: zh-qs Date: Thu, 19 Jan 2023 10:08:11 +0100 Subject: [PATCH 02/13] refactor performance tests --- engine/CMakeLists.txt | 4 +- engine/src/Performance.cu | 110 ++++++++++++++++++++ engine/src/Performance.cuh | 37 +++++++ engine/src/main.cu | 96 ++++++++++++++---- engine/test/AdvancedIntegral.cu | 173 +++++--------------------------- 5 files changed, 250 insertions(+), 170 deletions(-) create mode 100644 engine/src/Performance.cu create mode 100644 engine/src/Performance.cuh diff --git a/engine/CMakeLists.txt b/engine/CMakeLists.txt index 7e47911..72abac3 100644 --- a/engine/CMakeLists.txt +++ b/engine/CMakeLists.txt @@ -72,7 +72,9 @@ add_library( ${PARSER_DIR}/Scanner.cu ${PARSER_DIR}/Parser.cu) -add_executable(cusymint ${SRC_DIR}/main.cu) +add_executable(cusymint + ${SRC_DIR}/Performance.cu + ${SRC_DIR}/main.cu) add_executable( testcusymint diff --git a/engine/src/Performance.cu b/engine/src/Performance.cu new file mode 100644 index 0000000..7e7def6 --- /dev/null +++ b/engine/src/Performance.cu @@ -0,0 +1,110 @@ +#include "Performance.cuh" + +namespace Performance { + namespace { + std::string exec_and_read_output(const std::string cmd) { + // sciagniete z SolverProcessManager.cu + std::array buffer{}; + std::string result; + std::unique_ptr pipe(popen(cmd.c_str(), "r"), pclose); + if (!pipe) { + throw std::runtime_error("popen() failed!"); + } + while (fgets(buffer.data(), buffer.size(), pipe.get()) != nullptr) { + result += buffer.data(); + } + return result; + } + + void capitalize_substring(std::string& string, const std::string& what) { + if (string.size() < what.size()) { + return; + } + for (size_t i = 0; i <= string.size() - what.size(); ++i) { + if (string.substr(i, what.size()) == what && + (i == 0 || !std::isalpha(string[i - 1])) && + (i == string.size() - what.size() || !std::isalpha(string[i + what.size()]))) { + string[i] = string[i] - 'a' + 'A'; + } + } + } + + void replace(std::string& string, const std::string& from, const std::string& to) { + if (string.size() < from.size()) { + return; + } + for (ssize_t i = 0; i <= string.size() - from.size(); ++i) { + if (string.substr(i, from.size()) == from) { + string = string.substr(0, i) + to + string.substr(i + from.size()); + } + } + } + + std::string make_mathematica_command(const std::string& integral) { + // /media/cmonbrug/DATA/mathematica/Executables/wolframscript -code + // 'expr=ToExpression["sin(Pi)+x",TraditionalForm];{t,b}=AbsoluteTiming[Integrate[expr,x]];t' + std::string expression = integral; + capitalize_substring(expression, "e"); + capitalize_substring(expression, "pi"); + capitalize_substring(expression, "sqrt"); + + return fmt::format( + R"(wolframscript -code 'expr=ToExpression["{}",TraditionalForm];{{t,b}}=AbsoluteTiming[Integrate[expr,x]];t')", + expression); + } + + std::string make_sympy_command(const std::string& integral) { + std::string expression = integral; + capitalize_substring(expression, "e"); + replace(expression, "^", "**"); + replace(expression, "arctan", "atan"); + replace(expression, "arccos", "acos"); + replace(expression, "arcsin", "asin"); + replace(expression, "sgn", "sign"); + + return fmt::format( + R"(python3 -c 'from sympy import *;x=Symbol("x");t=Symbol("t");print(utilities.timeutils.timed(lambda:integrate({},x))[1])')", + expression); + } + } + + void do_not_print_results(const std::string& integral_str, const double& cusymint_seconds, + bool cusymint_success, const std::string& mathematica_result, + const std::string& sympy_result) {} + + void print_human_readable_results(const std::string& integral_str, + const double& cusymint_seconds, bool cusymint_success, + const std::string& mathematica_result, + const std::string& sympy_result) { + printf("%s:\n" + " Cusymint time: %f%s\n" + " Mathematica time: %s\n" + " Sympy time: %s\n", + integral_str.c_str(), cusymint_seconds, cusymint_success ? "" : " (failure)", + mathematica_result.c_str(), sympy_result.c_str()); + } + + void print_csv_results(const std::string& integral_str, const double& cusymint_seconds, + bool cusymint_success, const std::string& mathematica_result, + const std::string& sympy_result) { + printf("%s;%s;%f;%s;%s\n", integral_str.c_str(), cusymint_success ? "TRUE" : "FALSE", + cusymint_seconds, mathematica_result.c_str(), sympy_result.c_str()); + } + + void test_with_other_solutions(const std::string& integral_str, PrintRoutine print_results) { + const auto integral = Sym::integral(Parser::parse_function(integral_str)); + + Sym::Integrator integrator; + + const clock_t start = clock(); + const auto result = integrator.solve_integral(integral); + const clock_t end = clock(); + + const double cusymint_seconds = static_cast(end - start) / CLOCKS_PER_SEC; + + auto mathematica_result = exec_and_read_output(make_mathematica_command(integral_str)); + auto sympy_result = exec_and_read_output(make_sympy_command(integral[1].to_string())); + + print_results(integral_str, cusymint_seconds, result.has_value(), mathematica_result, sympy_result); + } +} \ No newline at end of file diff --git a/engine/src/Performance.cuh b/engine/src/Performance.cuh new file mode 100644 index 0000000..58c777d --- /dev/null +++ b/engine/src/Performance.cuh @@ -0,0 +1,37 @@ +#ifndef PERFORMACNE_CUH +#define PERFORMANCE_CUH + +#include +#include +#include +#include +#include + +#include "Evaluation/Integrator.cuh" +#include "Parser/Parser.cuh" +#include "Symbol/Macros.cuh" +#include "Symbol/Symbol.cuh" + +namespace Performance { + + using PrintRoutine = void (*)(const std::string& integral_str, const double& cusymint_seconds, + bool cusymint_success, const std::string& mathematica_result, + const std::string& sympy_result); + + void do_not_print_results(const std::string& integral_str, const double& cusymint_seconds, + bool cusymint_success, const std::string& mathematica_result, + const std::string& sympy_result); + void print_human_readable_results(const std::string& integral_str, + const double& cusymint_seconds, bool cusymint_success, + const std::string& mathematica_result, + const std::string& sympy_result); + + void print_csv_results(const std::string& integral_str, const double& cusymint_seconds, + bool cusymint_success, const std::string& mathematica_result, + const std::string& sympy_result); + + void test_with_other_solutions(const std::string& integral_str, + PrintRoutine print_results = print_human_readable_results); +} + +#endif \ No newline at end of file diff --git a/engine/src/main.cu b/engine/src/main.cu index d5b9a56..09c987e 100644 --- a/engine/src/main.cu +++ b/engine/src/main.cu @@ -1,6 +1,5 @@ #include - #include #include #include @@ -10,6 +9,7 @@ #include "Evaluation/Integrator.cuh" #include "Evaluation/StaticFunctions.cuh" +#include "Performance.cuh" #include "Symbol/Constants.cuh" #include "Symbol/ExpressionArray.cuh" #include "Symbol/Integral.cuh" @@ -43,6 +43,18 @@ std::string e_tower(size_t n) { return res; } +void cuda_warmup() { + Sym::Integrator integrator; + integrator.solve_integral(Sym::integral(Sym::var() + Sym::sin(Sym::var()))); +} + +void test_performance(const std::vector& integrals) { + for (auto integral_str : integrals) { + Performance::test_with_other_solutions(integral_str, + Performance::print_human_readable_results); + } +} + int main() { if constexpr (Consts::DEBUG) { fmt::print("Running in debug mode\n"); @@ -50,25 +62,65 @@ int main() { Sym::Static::init_functions(); - - const auto integral = Sym::integral(Parser::parse_function("tg(x)^2")); - - fmt::print("Trying to solve an integral: {}\n", integral.data()->to_tex()); - - Sym::Integrator integrator; - Sym::ComputationHistory history; - const auto solution = integrator.solve_integral_with_history(integral, history); - - if (solution.has_value()) { - fmt::print("Success! Solution:\n{} + C\n", solution.value().data()->to_tex()); - - fmt::print("\nComputation steps:\n\n"); - - for (const auto& step: history.get_tex_history()) { - fmt::print(" {}\\\\\n\n", step); - } - } - else { - fmt::print("No solution found\n"); - } + cuda_warmup(); + + test_performance({ + "x^4-5x+3", + "(1-x^3)^2", + "x^3t^2", + "x^(2/3)", + "1/x^(2/3)", + "(x^3-3x^2+1)/x^5", + "(x^(2/3)-sqrt(x))/x^(1/3)", + "2e^x+3*5^x", + "x cos(x)", + "x e^x", + "x^2e^x", + "x^2cos(x)", + "sqrt(x)ln(x)", + "ln(x)/x^4", + "(2x-1)^20", + "cos(3x-1)", + "e^x/(3+e^x)", + "tg(x)", + "tg(x)/cos^2(x)", + "cos(x)e^sin(x)", + "cos(x)/sqrt(1+sin(x))", + "6^(1-x)", + "(sqrt(x)-2)^2/x^2", + "1/(sin^2(x)cos^2(x))", + "x*arctg(x)", + "abs(x)", + "(abs(1-x^2)+1+x^2)/2", + "tg^2(x)", + "(2x-3)^10", + "sin^5(x)cos(x)", + "ln(x)", + "x*cos(x)", + "x^2e^(1-x)", + "sin^5(x)", + "sin^4(x)cos^3(x)", + "e^x*e^e^x*e^e^e^x*e^e^e^e^x*e^e^e^e^e^x*e^e^e^e^e^e^x", + }); + + // const auto integral = Sym::integral(Parser::parse_function("tg(x)^2")); + + // fmt::print("Trying to solve an integral: {}\n", integral.data()->to_tex()); + + // Sym::Integrator integrator; + // Sym::ComputationHistory history; + // const auto solution = integrator.solve_integral_with_history(integral, history); + + // if (solution.has_value()) { + // fmt::print("Success! Solution:\n{} + C\n", solution.value().data()->to_tex()); + + // fmt::print("\nComputation steps:\n\n"); + + // for (const auto& step: history.get_tex_history()) { + // fmt::print(" {}\\\\\n\n", step); + // } + // } + // else { + // fmt::print("No solution found\n"); + // } } diff --git a/engine/test/AdvancedIntegral.cu b/engine/test/AdvancedIntegral.cu index 0a4cc96..854f7cb 100644 --- a/engine/test/AdvancedIntegral.cu +++ b/engine/test/AdvancedIntegral.cu @@ -1,114 +1,11 @@ #include "IntegralCommons.cuh" -#include -#include -#include -#include -#include - #include "Evaluation/Integrator.cuh" -#include "Parser/Parser.cuh" -#include "Simplify.cuh" -#include "Symbol/Macros.cuh" -#include "Symbol/Symbol.cuh" - -#define ADVANCED_INTEGRAL_TEST_CUDA_ONLY(_name, _integral, _expected_result) \ - INTEGRAL_TEST(AdvancedIntegral, _name, _integral, _expected_result) - -#define ADVANCED_INTEGRAL_TEST_TIME(_name, _integral, _expected_result) \ - TEST_F(AdvancedIntegral, _name) { test_integral_with_external_tools(_integral); } #define ADVANCED_INTEGRAL_TEST(_name, _integral, _expected_result) \ - ADVANCED_INTEGRAL_TEST_TIME(_name, _integral, _expected_result) + INTEGRAL_TEST(AdvancedIntegral, _name, _integral, _expected_result) namespace Test { - namespace { - std::string exec_and_read_output(const std::string cmd) { - // sciagniete z SolverProcessManager.cu - std::array buffer{}; - std::string result; - std::unique_ptr pipe(popen(cmd.c_str(), "r"), pclose); - if (!pipe) { - throw std::runtime_error("popen() failed!"); - } - while (fgets(buffer.data(), buffer.size(), pipe.get()) != nullptr) { - result += buffer.data(); - } - return result; - } - - void capitalize_substring(std::string& string, const std::string& what) { - if (string.size() < what.size()) { - return; - } - for (size_t i = 0; i <= string.size() - what.size(); ++i) { - if (string.substr(i, what.size()) == what && - (i == 0 || !std::isalpha(string[i - 1])) && - (i == string.size() - what.size() || !std::isalpha(string[i + what.size()]))) { - string[i] = string[i] - 'a' + 'A'; - } - } - } - - void replace(std::string& string, const std::string& from, const std::string& to) { - if (string.size() < from.size()) { - return; - } - for (ssize_t i = 0; i <= string.size() - from.size(); ++i) { - if (string.substr(i, from.size()) == from) { - string = string.substr(0, i) + to + string.substr(i + from.size()); - } - } - } - - std::string make_mathematica_command(const std::string& integral) { - // /media/cmonbrug/DATA/mathematica/Executables/wolframscript -code - // 'expr=ToExpression["sin(Pi)+x",TraditionalForm];{t,b}=AbsoluteTiming[Integrate[expr,x]];t' - std::string expression = integral; - capitalize_substring(expression, "e"); - capitalize_substring(expression, "pi"); - capitalize_substring(expression, "sqrt"); - return fmt::format( - R"(wolframscript -code 'expr=ToExpression["{}",TraditionalForm];{{t,b}}=AbsoluteTiming[Integrate[expr,x]];t')", - expression); - } - - std::string make_sympy_command(const std::string& integral) { - std::string expression = integral; - capitalize_substring(expression, "e"); - replace(expression, "^", "**"); - replace(expression, "arctan", "atan"); - replace(expression, "arccos", "acos"); - replace(expression, "arcsin", "asin"); - replace(expression, "sgn", "sign"); - //fmt::print(R"(Call: python3 -c 'from sympy import *;x=Symbol("x");print(utilities.timeutils.timed(lambda:integrate({},x))[1])')" "\n", - // expression); - return fmt::format( - R"(python3 -c 'from sympy import *;x=Symbol("x");t=Symbol("t");print(utilities.timeutils.timed(lambda:integrate({},x))[1])')", - expression); - } - - void test_integral_with_external_tools(const std::string& integral_str) { - const auto integral = Sym::integral(Parser::parse_function(integral_str)); - - Sym::Integrator integrator; - - const clock_t start = clock(); - const auto result = integrator.solve_integral(integral); - const clock_t end = clock(); - - const double cusymint_seconds = static_cast(end - start) / CLOCKS_PER_SEC; - - auto mathematica_result = exec_and_read_output(make_mathematica_command(integral_str)); - auto sympy_result = exec_and_read_output(make_sympy_command(integral[1].to_string())); - - printf("Cusymint time: %f\n" - "Mathematica time: %s\n" - "Sympy time: %s\n", - cusymint_seconds, mathematica_result.c_str(), sympy_result.c_str()); - } - } - class AdvancedIntegral : public IntegrationFixture {}; // easy excercises from http://www.math.uni.wroc.pl/~ikrol/zbiorek7.pdf @@ -124,29 +21,24 @@ namespace Test { ADVANCED_INTEGRAL_TEST(CosineByParts, "x cos(x)", "x sin(x)+cos(x)") ADVANCED_INTEGRAL_TEST(ExponentByParts, "x e^x", "x e^x-e^x") ADVANCED_INTEGRAL_TEST(ExponentByParts2, "x^2e^x", "e^x(x^2-2x+2)") - // ADVANCED_INTEGRAL_TEST(NegativeExponentByParts, "x e^(-x)", "-e^(-x)(x+1)") // linear - // substitution for nested-vars-only required + ADVANCED_INTEGRAL_TEST(NegativeExponentByParts, "x e^(-x)", "-e^(-x)(x+1)") ADVANCED_INTEGRAL_TEST(CosineByParts2, "x^2cos(x)", "x^2sin(x)+2x cos(x)-2sin(x)") ADVANCED_INTEGRAL_TEST(LogarithmTimesRoot, "sqrt(x)ln(x)", "2/3x^(3/2)(ln(x)-2/3)") ADVANCED_INTEGRAL_TEST(LogarithmDividedByPower, "ln(x)/x^4", "-1/(3x^3)(ln(x)+1/3)") ADVANCED_INTEGRAL_TEST(VeryLongPolynomial, "(2x-1)^20", "1/42(2x-1)^21+1/42") ADVANCED_INTEGRAL_TEST(CosineOfLinearFunction, "cos(3x-1)", "sin(3x-1)/3") - // ADVANCED_INTEGRAL_TEST(RationalWithSubstitution, "2x/(5x^2+1)", "1/5ln(5x^2+1)") // - // substitution t=x^2+a required + // ADVANCED_INTEGRAL_TEST(RationalWithSubstitution, "2x/(5x^2+1)", "1/5ln(5x^2+1)") // substitution t=x^2+a required ADVANCED_INTEGRAL_TEST(SimpleRationalWithExponent, "e^x/(3+e^x)", "ln(3+e^x)") - // ADVANCED_INTEGRAL_TEST(ExponentWithSquare, "x*e^(-x^2)", "-1/2e^(-x^2)") // substitution - // t=x^2+a required ADVANCED_INTEGRAL_TEST(ExponentOfInverse, "e^(1/x)/x^2", "-e^(1/x)") // - // substitution t=ax^b+c required (dx = 1/b((t-c)/a)^((1-b)/b)) + // ADVANCED_INTEGRAL_TEST(ExponentWithSquare, "x*e^(-x^2)", "-1/2e^(-x^2)") // substitution t=x^2+a required + // ADVANCED_INTEGRAL_TEST(ExponentOfInverse, "e^(1/x)/x^2", "-e^(1/x)") // substitution t=ax^b+c required (dx = 1/b((t-c)/a)^((1-b)/b)) ADVANCED_INTEGRAL_TEST(Tangent, "tg(x)", "-ln(cos(x))") ADVANCED_INTEGRAL_TEST(TangentDividedByCosSquared, "tg(x)/cos^2(x)", "tg^2(x)/2") ADVANCED_INTEGRAL_TEST(CosTimesExpSine, "cos(x)e^sin(x)", "e^sin(x)") ADVANCED_INTEGRAL_TEST(CosineDividedByRootOfLinSine, "cos(x)/sqrt(1+sin(x))", "2sqrt(1+sin(x))") - // ADVANCED_INTEGRAL_TEST(TangentDerivativeWithSubstitution, "x^3/cos^2(x^4)", "tg(x^4)/4") // - // substitution t=ax^b+c required (dx = 1/b((t-c)/a)^((1-b)/b)) + // ADVANCED_INTEGRAL_TEST(TangentDerivativeWithSubstitution, "x^3/cos^2(x^4)", "tg(x^4)/4") // substitution t=ax^b+c required (dx = 1/b((t-c)/a)^((1-b)/b)) ADVANCED_INTEGRAL_TEST(ExponentOfLinear, "6^(1-x)", "-6^(1-x)/ln(6)") - // ADVANCED_INTEGRAL_TEST(ArctanSubtitution, "1/(1+x^2)/arctan(x)", "ln(arctan(x))") // - // substitution t=arctan(x) required + // ADVANCED_INTEGRAL_TEST(ArctanSubtitution, "1/(1+x^2)/arctan(x)", "ln(arctan(x))") // substitution t=arctan(x) required // task 4 contains examples which require nested integrals (not implemeneted) @@ -155,21 +47,16 @@ namespace Test { ADVANCED_INTEGRAL_TEST(SimpleRationalWithSquare, "(sqrt(x)-2)^2/x^2", "-4/x+8/sqrt(x)+ln(x)") ADVANCED_INTEGRAL_TEST(SineCosineSquaredInDenominator, "1/(sin^2(x)cos^2(x))", "-1/tg(x)+tg(x)") - // ADVANCED_INTEGRAL_TEST(LogarithmDividedByX, "ln^5(x)/x", "ln^6(x)/6") // substitution t=ln(x) - // or cyclic integrals required ADVANCED_INTEGRAL_TEST(ArctangentWithX3Substitution, - // "x^3/(x^8+1)", "arctg(x^4)/4") // substitution t=x^4 required - // ADVANCED_INTEGRAL_TEST(LogarithmMultipliedByX, "x^n*ln(x)", - // "x^(n+1)ln(x)/(n+1)-x^(n+1)/(n+1)^2") // non-numeric powers integration by parts required + // ADVANCED_INTEGRAL_TEST(LogarithmDividedByX, "ln^5(x)/x", "ln^6(x)/6") // substitution t=ln(x) or cyclic integrals required + // ADVANCED_INTEGRAL_TEST(ArctangentWithX3Substitution, "x^3/(x^8+1)", "arctg(x^4)/4") // substitution t=x^4 required + // ADVANCED_INTEGRAL_TEST(LogarithmMultipliedByX, "x^n*ln(x)", "x^(n+1)ln(x)/(n+1)-x^(n+1)/(n+1)^2") // non-numeric powers integration by parts required ADVANCED_INTEGRAL_TEST(ArctangentMultipliedByX, "x*arctg(x)", "x^2/2*arctg(x)-x/2+arctg(x)/2") - // ADVANCED_INTEGRAL_TEST(Arcsine, "arcsin(x)", "x*arcsin(x)+2*sqrt(1-x^2)") // substitution - // t=x^2+a required ADVANCED_INTEGRAL_TEST(LogarithmWithX2Substitution, "x*ln(x^2+1)", - // "(x^2+1)ln(x^2+1)-x^2-1") // substitution t=x^2+a required - // ADVANCED_INTEGRAL_TEST(CyclicIntegral1, "e^x*sin(2x)", "e^x*sin(2x)/5 - e^x*2cos(2x)/5") // - // cyclic integrals required ADVANCED_INTEGRAL_TEST(CyclicIntegral2, "sin(ln(x))", - // "-x/2(cos(ln(x))-sin(ln(x)))") // cyclic integrals required - // ADVANCED_INTEGRAL_TEST(CotangentDividedByLogOfSine, "ctg(x)/ln(sin(x))", "ln(ln(sin(x)))") // - // substitution t=ln(x) required ADVANCED_INTEGRAL_TEST(InvertedCosH, "1/(e^x+e^-x)", - // "arctg(e^x)") // extract_function required + // ADVANCED_INTEGRAL_TEST(Arcsine, "arcsin(x)", "x*arcsin(x)+2*sqrt(1-x^2)") // substitution t=x^2+a required + // ADVANCED_INTEGRAL_TEST(LogarithmWithX2Substitution, "x*ln(x^2+1)", "(x^2+1)ln(x^2+1)-x^2-1") // substitution t=x^2+a required + // ADVANCED_INTEGRAL_TEST(CyclicIntegral1, "e^x*sin(2x)", "e^x*sin(2x)/5 - e^x*2cos(2x)/5") // cyclic integrals required + // ADVANCED_INTEGRAL_TEST(CyclicIntegral2, "sin(ln(x))", "-x/2(cos(ln(x))-sin(ln(x)))") // cyclic integrals required + // ADVANCED_INTEGRAL_TEST(CotangentDividedByLogOfSine, "ctg(x)/ln(sin(x))", "ln(ln(sin(x)))") // substitution t=ln(x) required + // ADVANCED_INTEGRAL_TEST(InvertedCosH, "1/(e^x+e^-x)", "arctg(e^x)") // extract_function required ADVANCED_INTEGRAL_TEST(Absolute, "abs(x)", "x*abs(x)/2") // abs required ADVANCED_INTEGRAL_TEST(MaxOfOneAndSquare, "(abs(1-x^2)+1+x^2)/2", @@ -180,30 +67,22 @@ namespace Test { // https://pages.mini.pw.edu.pl/~dembinskaa/www/?download=Inf_I_PowtKol3_2022-2023.pdf, task 11 ADVANCED_INTEGRAL_TEST(TangentSquared, "tg^2(x)", "tg(x)-arctg(tg(x))") // may sometimes fail ADVANCED_INTEGRAL_TEST(LongPolynomial, "(2x-3)^10", "1/22(2x-3)^11+3^11/22") - // ADVANCED_INTEGRAL_TEST(SquareRootInDenominator, "1/(2+sqrt(x))", "2sqrt(x)-4ln(sqrt(x)+2)") - // // substitution t=sqrt(x)+a required ADVANCED_INTEGRAL_TEST(InvertedCosH2, "1/cosh(x)", - // "2arctg(e^x)") // extract_funcion required + // ADVANCED_INTEGRAL_TEST(SquareRootInDenominator, "1/(2+sqrt(x))", "2sqrt(x)-4ln(sqrt(x)+2)") // substitution t=sqrt(x)+a required + // ADVANCED_INTEGRAL_TEST(InvertedCosH2, "1/cosh(x)", "2arctg(e^x)") // extract_funcion required ADVANCED_INTEGRAL_TEST(Sine5Cosine, "sin^5(x)cos(x)", "sin^6(x)/6") - // ADVANCED_INTEGRAL_TEST(ExpressionWithSquareInDenominator, "x/(x^2-1)^(3/2)", - // "-1/sqrt(x^2-1)") // substitution t=ax^b+c required - // ADVANCED_INTEGRAL_TEST(ArcsineWithDerivative, "arcsin^2(x)/sqrt(1-x^2)", "arcsin^3(x)/3") // - // substitution t=arcsin(x) required ADVANCED_INTEGRAL_TEST(RootInDenominator, "1/(x^(1/3)+1)", - // "3/2*x^(2/3)-3x^(1/3)+3ln(x^(1/3)+1)") // substitution t=ax^b+c required (dx = - // 1/b((t-c)/a)^((1-b)/b)) + // ADVANCED_INTEGRAL_TEST(ExpressionWithSquareInDenominator, "x/(x^2-1)^(3/2)", "-1/sqrt(x^2-1)") // substitution t=ax^b+c required + // ADVANCED_INTEGRAL_TEST(ArcsineWithDerivative, "arcsin^2(x)/sqrt(1-x^2)", "arcsin^3(x)/3") // substitution t=arcsin(x) required + // ADVANCED_INTEGRAL_TEST(RootInDenominator, "1/(x^(1/3)+1)", "3/2*x^(2/3)-3x^(1/3)+3ln(x^(1/3)+1)") // substitution t=ax^b+c required (dx = 1/b((t-c)/a)^((1-b)/b)) ADVANCED_INTEGRAL_TEST(Logarithm, "ln(x)", "x ln(x) - x") ADVANCED_INTEGRAL_TEST(XTimesCosine, "x*cos(x)", "x*sin(x)+cos(x)") - // ADVANCED_INTEGRAL_TEST(XTimesExponential, "x^2e^(1-x)", "-(x^2+2x+2)e^(1-x)") // linear - // substitution for nested-vars-only required ADVANCED_INTEGRAL_TEST(SquareTimesExponential, - // "8x^2e^(4-x^3)", "-8/3e^(4-x^3)") // substitution t=x^2+a required + ADVANCED_INTEGRAL_TEST(XTimesExponential, "x^2e^(1-x)", "-(x^2+2x+2)e^(1-x)") + // ADVANCED_INTEGRAL_TEST(SquareTimesExponential, "8x^2e^(4-x^3)", "-8/3e^(4-x^3)") // substitution t=x^2+a required - // ADVANCED_INTEGRAL_TEST(SineSquared, "sin^2(x)", "x/2-sin(2x)/4") // rational integrals - // required + // ADVANCED_INTEGRAL_TEST(SineSquared, "sin^2(x)", "x/2-sin(2x)/4") // rational integrals required ADVANCED_INTEGRAL_TEST(Sine5, "sin^5(x)", "-cos^5(x)/5+2/3cos^3(x)-cos(x)") ADVANCED_INTEGRAL_TEST(Sine4Cos3, "sin^4(x)cos^3(x)", "sin^5(x)/5-sin^7(x)/7") - // ADVANCED_INTEGRAL_TEST(Sine4Cos2, "sin^4(x)cos^2(x)", - // "1/6sin^5(x)cos(x)+1/24sin^3(x)cos(x)-1/32sin(x)cos(x)+1/16x") // rational integrals required - // ADVANCED_INTEGRAL_TEST(Cos3XCos5X, "cos(3x)cos(5x)", "1/16sin(8x)+1/4sin(2x)") // more - // trigonometric identitied required + // ADVANCED_INTEGRAL_TEST(Sine4Cos2, "sin^4(x)cos^2(x)", "1/6sin^5(x)cos(x)+1/24sin^3(x)cos(x)-1/32sin(x)cos(x)+1/16x") // rational integrals required + // ADVANCED_INTEGRAL_TEST(Cos3XCos5X, "cos(3x)cos(5x)", "1/16sin(8x)+1/4sin(2x)") // more trigonometric identitied required // own tests ADVANCED_INTEGRAL_TEST(ETower, "e^x*e^e^x*e^e^e^x*e^e^e^e^x*e^e^e^e^e^x*e^e^e^e^e^e^x", From c342ce7b04dac3efcde3c6d6a71225f0075fb0f0 Mon Sep 17 00:00:00 2001 From: zh-qs Date: Thu, 19 Jan 2023 10:11:15 +0100 Subject: [PATCH 03/13] add +C in the last step --- engine/src/Evaluation/ComputationHistory.cu | 13 ++- engine/src/main.cu | 118 ++++++++++---------- 2 files changed, 68 insertions(+), 63 deletions(-) diff --git a/engine/src/Evaluation/ComputationHistory.cu b/engine/src/Evaluation/ComputationHistory.cu index e53ac6c..9118ed4 100644 --- a/engine/src/Evaluation/ComputationHistory.cu +++ b/engine/src/Evaluation/ComputationHistory.cu @@ -290,19 +290,24 @@ namespace Sym { else { size_t trans_idx = 1; for (const auto& trans : step.get_operations(*prev_step)) { - result.push_back(fmt::format( - R"(\quad \text{{{}{} {}:}}\:{}:)", trans_idx, get_ordinal_suffix(trans_idx), - Globalization::INTEGRAL, trans->get_description())); + result.push_back(fmt::format(R"(\quad \text{{{}{} {}:}}\:{}:)", trans_idx, + get_ordinal_suffix(trans_idx), + Globalization::INTEGRAL, + trans->get_description())); ++trans_idx; } } - result.push_back(fmt::format(R"(=\qquad {})", expression.data()->to_tex())); + result.push_back(fmt::format(R"(=\qquad {})", expression_str)); } prev_str = expression_str; prev_step = &step; } + if (!result.empty()) { + result[result.size() - 1] += " + C"; + } + return result; } diff --git a/engine/src/main.cu b/engine/src/main.cu index 09c987e..312ee56 100644 --- a/engine/src/main.cu +++ b/engine/src/main.cu @@ -64,63 +64,63 @@ int main() { cuda_warmup(); - test_performance({ - "x^4-5x+3", - "(1-x^3)^2", - "x^3t^2", - "x^(2/3)", - "1/x^(2/3)", - "(x^3-3x^2+1)/x^5", - "(x^(2/3)-sqrt(x))/x^(1/3)", - "2e^x+3*5^x", - "x cos(x)", - "x e^x", - "x^2e^x", - "x^2cos(x)", - "sqrt(x)ln(x)", - "ln(x)/x^4", - "(2x-1)^20", - "cos(3x-1)", - "e^x/(3+e^x)", - "tg(x)", - "tg(x)/cos^2(x)", - "cos(x)e^sin(x)", - "cos(x)/sqrt(1+sin(x))", - "6^(1-x)", - "(sqrt(x)-2)^2/x^2", - "1/(sin^2(x)cos^2(x))", - "x*arctg(x)", - "abs(x)", - "(abs(1-x^2)+1+x^2)/2", - "tg^2(x)", - "(2x-3)^10", - "sin^5(x)cos(x)", - "ln(x)", - "x*cos(x)", - "x^2e^(1-x)", - "sin^5(x)", - "sin^4(x)cos^3(x)", - "e^x*e^e^x*e^e^e^x*e^e^e^e^x*e^e^e^e^e^x*e^e^e^e^e^e^x", - }); - - // const auto integral = Sym::integral(Parser::parse_function("tg(x)^2")); - - // fmt::print("Trying to solve an integral: {}\n", integral.data()->to_tex()); - - // Sym::Integrator integrator; - // Sym::ComputationHistory history; - // const auto solution = integrator.solve_integral_with_history(integral, history); - - // if (solution.has_value()) { - // fmt::print("Success! Solution:\n{} + C\n", solution.value().data()->to_tex()); - - // fmt::print("\nComputation steps:\n\n"); - - // for (const auto& step: history.get_tex_history()) { - // fmt::print(" {}\\\\\n\n", step); - // } - // } - // else { - // fmt::print("No solution found\n"); - // } + // test_performance({ + // "x^4-5x+3", + // "(1-x^3)^2", + // "x^3t^2", + // "x^(2/3)", + // "1/x^(2/3)", + // "(x^3-3x^2+1)/x^5", + // "(x^(2/3)-sqrt(x))/x^(1/3)", + // "2e^x+3*5^x", + // "x cos(x)", + // "x e^x", + // "x^2e^x", + // "x^2cos(x)", + // "sqrt(x)ln(x)", + // "ln(x)/x^4", + // "(2x-1)^20", + // "cos(3x-1)", + // "e^x/(3+e^x)", + // "tg(x)", + // "tg(x)/cos^2(x)", + // "cos(x)e^sin(x)", + // "cos(x)/sqrt(1+sin(x))", + // "6^(1-x)", + // "(sqrt(x)-2)^2/x^2", + // "1/(sin^2(x)cos^2(x))", + // "x*arctg(x)", + // "abs(x)", + // "(abs(1-x^2)+1+x^2)/2", + // "tg^2(x)", + // "(2x-3)^10", + // "sin^5(x)cos(x)", + // "ln(x)", + // "x*cos(x)", + // "x^2e^(1-x)", + // "sin^5(x)", + // "sin^4(x)cos^3(x)", + // "e^x*e^e^x*e^e^e^x*e^e^e^e^x*e^e^e^e^e^x*e^e^e^e^e^e^x", + // }); + + const auto integral = Sym::integral(Parser::parse_function("tg(x)^2")); + + fmt::print("Trying to solve an integral: {}\n", integral.data()->to_tex()); + + Sym::Integrator integrator; + Sym::ComputationHistory history; + const auto solution = integrator.solve_integral_with_history(integral, history); + + if (solution.has_value()) { + fmt::print("Success! Solution:\n{} + C\n", solution.value().data()->to_tex()); + + fmt::print("\nComputation steps:\n\n"); + + for (const auto& step: history.get_tex_history()) { + fmt::print(" {}\\\\\n\n", step); + } + } + else { + fmt::print("No solution found\n"); + } } From 553b6bfa7bab500f4c21cbb951fb65c27201cbe8 Mon Sep 17 00:00:00 2001 From: zh-qs Date: Thu, 19 Jan 2023 11:06:54 +0100 Subject: [PATCH 04/13] add matlab --- engine/src/Performance.cu | 38 +++++++---- engine/src/Performance.cuh | 8 +-- engine/src/main.cu | 128 ++++++++++++++++++------------------- 3 files changed, 95 insertions(+), 79 deletions(-) diff --git a/engine/src/Performance.cu b/engine/src/Performance.cu index 7e7def6..ce1be43 100644 --- a/engine/src/Performance.cu +++ b/engine/src/Performance.cu @@ -49,7 +49,7 @@ namespace Performance { capitalize_substring(expression, "sqrt"); return fmt::format( - R"(wolframscript -code 'expr=ToExpression["{}",TraditionalForm];{{t,b}}=AbsoluteTiming[Integrate[expr,x]];t')", + R"(wolframscript -code 'expr=ToExpression["{}",TraditionalForm];Integrate[expr,x];{{t,b}}=AbsoluteTiming[Integrate[expr,x]];t')", expression); } @@ -63,32 +63,47 @@ namespace Performance { replace(expression, "sgn", "sign"); return fmt::format( - R"(python3 -c 'from sympy import *;x=Symbol("x");t=Symbol("t");print(utilities.timeutils.timed(lambda:integrate({},x))[1])')", + R"(python3 -c 'from sympy import *;x=Symbol("x");t=Symbol("t");integrate({},x);print(utilities.timeutils.timed(lambda:integrate({},x))[1])')", + expression, expression); + } + + std::string make_matlab_command(const std::string& integral) { + std::string expression = integral; + capitalize_substring(expression, "e"); + replace(expression, "E", "exp(sym(1))"); + replace(expression, "arctan", "atan"); + replace(expression, "arccos", "acos"); + replace(expression, "arcsin", "asin"); + replace(expression, "sgn", "sign"); + + return fmt::format( + R"(matlab -batch 'syms x z;f=@()int({},x);f();fprintf("%f\n",timeit(f));')", expression); } } void do_not_print_results(const std::string& integral_str, const double& cusymint_seconds, bool cusymint_success, const std::string& mathematica_result, - const std::string& sympy_result) {} + const std::string& sympy_result, const std::string& matlab_result) {} void print_human_readable_results(const std::string& integral_str, const double& cusymint_seconds, bool cusymint_success, const std::string& mathematica_result, - const std::string& sympy_result) { + const std::string& sympy_result, const std::string& matlab_result) { printf("%s:\n" " Cusymint time: %f%s\n" - " Mathematica time: %s\n" - " Sympy time: %s\n", + " Mathematica time: %s" + " Sympy time: %s" + " Matlab time: %s\n", integral_str.c_str(), cusymint_seconds, cusymint_success ? "" : " (failure)", - mathematica_result.c_str(), sympy_result.c_str()); + mathematica_result.c_str(), sympy_result.c_str(), matlab_result.c_str()); } void print_csv_results(const std::string& integral_str, const double& cusymint_seconds, bool cusymint_success, const std::string& mathematica_result, - const std::string& sympy_result) { - printf("%s;%s;%f;%s;%s\n", integral_str.c_str(), cusymint_success ? "TRUE" : "FALSE", - cusymint_seconds, mathematica_result.c_str(), sympy_result.c_str()); + const std::string& sympy_result, const std::string& matlab_result) { + printf("%s;%s;%f;%s;%s;%s\n", integral_str.c_str(), cusymint_success ? "TRUE" : "FALSE", + cusymint_seconds, mathematica_result.c_str(), sympy_result.c_str(), matlab_result.c_str()); } void test_with_other_solutions(const std::string& integral_str, PrintRoutine print_results) { @@ -104,7 +119,8 @@ namespace Performance { auto mathematica_result = exec_and_read_output(make_mathematica_command(integral_str)); auto sympy_result = exec_and_read_output(make_sympy_command(integral[1].to_string())); + auto matlab_result = exec_and_read_output(make_matlab_command(integral[1].to_string())); - print_results(integral_str, cusymint_seconds, result.has_value(), mathematica_result, sympy_result); + print_results(integral_str, cusymint_seconds, result.has_value(), mathematica_result, sympy_result, matlab_result); } } \ No newline at end of file diff --git a/engine/src/Performance.cuh b/engine/src/Performance.cuh index 58c777d..73786c2 100644 --- a/engine/src/Performance.cuh +++ b/engine/src/Performance.cuh @@ -16,19 +16,19 @@ namespace Performance { using PrintRoutine = void (*)(const std::string& integral_str, const double& cusymint_seconds, bool cusymint_success, const std::string& mathematica_result, - const std::string& sympy_result); + const std::string& sympy_result, const std::string& matlab_result); void do_not_print_results(const std::string& integral_str, const double& cusymint_seconds, bool cusymint_success, const std::string& mathematica_result, - const std::string& sympy_result); + const std::string& sympy_result, const std::string& matlab_result); void print_human_readable_results(const std::string& integral_str, const double& cusymint_seconds, bool cusymint_success, const std::string& mathematica_result, - const std::string& sympy_result); + const std::string& sympy_result, const std::string& matlab_result); void print_csv_results(const std::string& integral_str, const double& cusymint_seconds, bool cusymint_success, const std::string& mathematica_result, - const std::string& sympy_result); + const std::string& sympy_result, const std::string& matlab_result); void test_with_other_solutions(const std::string& integral_str, PrintRoutine print_results = print_human_readable_results); diff --git a/engine/src/main.cu b/engine/src/main.cu index 312ee56..e7dc672 100644 --- a/engine/src/main.cu +++ b/engine/src/main.cu @@ -43,9 +43,8 @@ std::string e_tower(size_t n) { return res; } -void cuda_warmup() { - Sym::Integrator integrator; - integrator.solve_integral(Sym::integral(Sym::var() + Sym::sin(Sym::var()))); +void warmup() { + Performance::test_with_other_solutions("x+sin(e^x)*e^x", Performance::do_not_print_results); } void test_performance(const std::vector& integrals) { @@ -62,65 +61,66 @@ int main() { Sym::Static::init_functions(); - cuda_warmup(); - - // test_performance({ - // "x^4-5x+3", - // "(1-x^3)^2", - // "x^3t^2", - // "x^(2/3)", - // "1/x^(2/3)", - // "(x^3-3x^2+1)/x^5", - // "(x^(2/3)-sqrt(x))/x^(1/3)", - // "2e^x+3*5^x", - // "x cos(x)", - // "x e^x", - // "x^2e^x", - // "x^2cos(x)", - // "sqrt(x)ln(x)", - // "ln(x)/x^4", - // "(2x-1)^20", - // "cos(3x-1)", - // "e^x/(3+e^x)", - // "tg(x)", - // "tg(x)/cos^2(x)", - // "cos(x)e^sin(x)", - // "cos(x)/sqrt(1+sin(x))", - // "6^(1-x)", - // "(sqrt(x)-2)^2/x^2", - // "1/(sin^2(x)cos^2(x))", - // "x*arctg(x)", - // "abs(x)", - // "(abs(1-x^2)+1+x^2)/2", - // "tg^2(x)", - // "(2x-3)^10", - // "sin^5(x)cos(x)", - // "ln(x)", - // "x*cos(x)", - // "x^2e^(1-x)", - // "sin^5(x)", - // "sin^4(x)cos^3(x)", - // "e^x*e^e^x*e^e^e^x*e^e^e^e^x*e^e^e^e^e^x*e^e^e^e^e^e^x", - // }); - - const auto integral = Sym::integral(Parser::parse_function("tg(x)^2")); - - fmt::print("Trying to solve an integral: {}\n", integral.data()->to_tex()); - - Sym::Integrator integrator; - Sym::ComputationHistory history; - const auto solution = integrator.solve_integral_with_history(integral, history); - - if (solution.has_value()) { - fmt::print("Success! Solution:\n{} + C\n", solution.value().data()->to_tex()); - - fmt::print("\nComputation steps:\n\n"); - - for (const auto& step: history.get_tex_history()) { - fmt::print(" {}\\\\\n\n", step); - } - } - else { - fmt::print("No solution found\n"); - } + printf("Warmup...\n\n"); + warmup(); + + test_performance({ + "x^4-5x+3", + "(1-x^3)^2", + "x^3t^2", + "x^(2/3)", + "1/x^(2/3)", + "(x^3-3x^2+1)/x^5", + "(x^(2/3)-sqrt(x))/x^(1/3)", + "2e^x+3*5^x", + "x cos(x)", + "x e^x", + "x^2e^x", + "x^2cos(x)", + "sqrt(x)ln(x)", + "ln(x)/x^4", + "(2x-1)^20", + "cos(3x-1)", + "e^x/(3+e^x)", + "tg(x)", + "tg(x)/cos^2(x)", + "cos(x)e^sin(x)", + "cos(x)/sqrt(1+sin(x))", + "6^(1-x)", + "(sqrt(x)-2)^2/x^2", + "1/(sin^2(x)cos^2(x))", + "x*arctg(x)", + "abs(x)", + "(abs(1-x^2)+1+x^2)/2", + "tg^2(x)", + "(2x-3)^10", + "sin^5(x)cos(x)", + "ln(x)", + "x*cos(x)", + "x^2e^(1-x)", + "sin^5(x)", + "sin^4(x)cos^3(x)", + "e^x*e^e^x*e^e^e^x*e^e^e^e^x*e^e^e^e^e^x*e^e^e^e^e^e^x", + }); + + // const auto integral = Sym::integral(Parser::parse_function("tg(x)^2")); + + // fmt::print("Trying to solve an integral: {}\n", integral.data()->to_tex()); + + // Sym::Integrator integrator; + // Sym::ComputationHistory history; + // const auto solution = integrator.solve_integral_with_history(integral, history); + + // if (solution.has_value()) { + // fmt::print("Success! Solution:\n{} + C\n", solution.value().data()->to_tex()); + + // fmt::print("\nComputation steps:\n\n"); + + // for (const auto& step: history.get_tex_history()) { + // fmt::print(" {}\\\\\n\n", step); + // } + // } + // else { + // fmt::print("No solution found\n"); + // } } From 681c0cf1aae067ceae388fa0a192e06510e4264f Mon Sep 17 00:00:00 2001 From: zh-qs Date: Thu, 19 Jan 2023 12:33:02 +0100 Subject: [PATCH 05/13] Make performance tests faster and more reliable --- engine/src/Performance.cu | 131 ++++++++++++++++++++++++++++++++++--- engine/src/Performance.cuh | 2 + engine/src/main.cu | 63 ++++++++---------- 3 files changed, 151 insertions(+), 45 deletions(-) diff --git a/engine/src/Performance.cu b/engine/src/Performance.cu index ce1be43..064ee05 100644 --- a/engine/src/Performance.cu +++ b/engine/src/Performance.cu @@ -1,9 +1,13 @@ #include "Performance.cuh" +#include "Symbol/Integral.cuh" +#include +#include +#include namespace Performance { namespace { std::string exec_and_read_output(const std::string cmd) { - // sciagniete z SolverProcessManager.cu + // stolen from SolverProcessManager.cu std::array buffer{}; std::string result; std::unique_ptr pipe(popen(cmd.c_str(), "r"), pclose); @@ -41,8 +45,6 @@ namespace Performance { } std::string make_mathematica_command(const std::string& integral) { - // /media/cmonbrug/DATA/mathematica/Executables/wolframscript -code - // 'expr=ToExpression["sin(Pi)+x",TraditionalForm];{t,b}=AbsoluteTiming[Integrate[expr,x]];t' std::string expression = integral; capitalize_substring(expression, "e"); capitalize_substring(expression, "pi"); @@ -80,6 +82,56 @@ namespace Performance { R"(matlab -batch 'syms x z;f=@()int({},x);f();fprintf("%f\n",timeit(f));')", expression); } + + std::string make_mathematica_command_batch(const std::vector& integrals) { + std::string result = "wolframscript -code \'"; + for (const auto& integral : integrals) { + std::string expression = integral; + capitalize_substring(expression, "e"); + capitalize_substring(expression, "pi"); + capitalize_substring(expression, "sqrt"); + result += fmt::format("{{t,b}}=AbsoluteTiming[Integrate[ToExpression[\"{}\"," + "TraditionalForm],x]];Print[t];", + expression); + } + + return result + "\'"; + } + + std::string make_sympy_command_batch(const std::vector& integrals) { + std::string result = R"(python3 -c 'from sympy import *;x=Symbol("x");t=Symbol("t");)"; + for (const auto& integral : integrals) { + std::string expression = integral; + capitalize_substring(expression, "e"); + replace(expression, "^", "**"); + replace(expression, "arctan", "atan"); + replace(expression, "arccos", "acos"); + replace(expression, "arcsin", "asin"); + replace(expression, "sgn", "sign"); + result += fmt::format( + "print(utilities.timeutils.timed(lambda:integrate({},x))[1]);", expression); + } + + return result + "\'"; + } + + std::string make_matlab_command_batch(const std::vector& integrals) { + std::string result = R"(matlab -batch 'syms x t;)"; + for (const auto& integral : integrals) { + std::string expression = integral; + capitalize_substring(expression, "e"); + replace(expression, "E", "exp(sym(1))"); + replace(expression, "arctan", "atan"); + replace(expression, "arccos", "acos"); + replace(expression, "arcsin", "asin"); + replace(expression, "sgn", "sign"); + replace(expression, "ln", "log"); + + result += fmt::format(R"(f=@()int({},x);fprintf("%f\n",timeit(f));)", expression); + } + + return result + "\'"; + } } void do_not_print_results(const std::string& integral_str, const double& cusymint_seconds, @@ -89,12 +141,13 @@ namespace Performance { void print_human_readable_results(const std::string& integral_str, const double& cusymint_seconds, bool cusymint_success, const std::string& mathematica_result, - const std::string& sympy_result, const std::string& matlab_result) { + const std::string& sympy_result, + const std::string& matlab_result) { printf("%s:\n" " Cusymint time: %f%s\n" - " Mathematica time: %s" - " Sympy time: %s" - " Matlab time: %s\n", + " Mathematica time: %s\n" + " Sympy time: %s\n" + " Matlab time: %s\n\n", integral_str.c_str(), cusymint_seconds, cusymint_success ? "" : " (failure)", mathematica_result.c_str(), sympy_result.c_str(), matlab_result.c_str()); } @@ -103,7 +156,8 @@ namespace Performance { bool cusymint_success, const std::string& mathematica_result, const std::string& sympy_result, const std::string& matlab_result) { printf("%s;%s;%f;%s;%s;%s\n", integral_str.c_str(), cusymint_success ? "TRUE" : "FALSE", - cusymint_seconds, mathematica_result.c_str(), sympy_result.c_str(), matlab_result.c_str()); + cusymint_seconds, mathematica_result.c_str(), sympy_result.c_str(), + matlab_result.c_str()); } void test_with_other_solutions(const std::string& integral_str, PrintRoutine print_results) { @@ -121,6 +175,65 @@ namespace Performance { auto sympy_result = exec_and_read_output(make_sympy_command(integral[1].to_string())); auto matlab_result = exec_and_read_output(make_matlab_command(integral[1].to_string())); - print_results(integral_str, cusymint_seconds, result.has_value(), mathematica_result, sympy_result, matlab_result); + print_results(integral_str, cusymint_seconds, result.has_value(), mathematica_result, + sympy_result, matlab_result); + } + + void test_performance(const std::vector& integrals, PrintRoutine print_results) { + std::vector cusymint_seconds_vector(integrals.size()); + std::vector cusymint_results_vector(integrals.size()); + std::vector mathematica_results_vector(integrals.size()); + std::vector sympy_results_vector(integrals.size()); + std::vector matlab_results_vector(integrals.size()); + + printf("Computing on CUDA...\n"); + + for (int i = 0; i < integrals.size(); ++i) { + Sym::Integrator integrator; + const clock_t start = clock(); + const auto result = + integrator.solve_integral(Sym::integral(Parser::parse_function(integrals[i]))); + const clock_t end = clock(); + + cusymint_seconds_vector[i] = static_cast(end - start) / CLOCKS_PER_SEC; + cusymint_results_vector[i] = result.has_value(); + } + + printf("Computing on Mathematica...\n"); + auto mathematica_result = exec_and_read_output(make_mathematica_command_batch(integrals)); + printf("Computing on SymPy...\n"); + auto sympy_result = exec_and_read_output(make_sympy_command_batch(integrals)); + printf("Computing on MATLAB...\n"); + auto matlab_result = exec_and_read_output(make_matlab_command_batch(integrals)); + + size_t mc_idx = 0; + size_t sp_idx = 0; + size_t ml_idx = 0; + for (int i = 0; i < integrals.size(); ++i) { + const auto mc_inc = mathematica_result.find('\n', mc_idx); + if (mc_inc != std::string::npos) { + mathematica_results_vector[i] = mathematica_result.substr(mc_idx, mc_inc - mc_idx); + mc_idx = mc_inc + 1; + } + + const auto sp_inc = sympy_result.find('\n', sp_idx); + if (sp_inc != std::string::npos) { + sympy_results_vector[i] = sympy_result.substr(sp_idx, sp_inc - sp_idx); + sp_idx = sp_inc + 1; + } + + const auto ml_inc = matlab_result.find('\n', ml_idx); + if (ml_inc != std::string::npos) { + matlab_results_vector[i] = matlab_result.substr(ml_idx, ml_inc - ml_idx); + ml_idx = ml_inc + 1; + } + } + + for (int i = 0; i < integrals.size(); ++i) { + print_results(integrals[i], cusymint_seconds_vector[i], cusymint_results_vector[i], + mathematica_results_vector[i], sympy_results_vector[i], + matlab_results_vector[i]); + } } + } \ No newline at end of file diff --git a/engine/src/Performance.cuh b/engine/src/Performance.cuh index 73786c2..d22fcba 100644 --- a/engine/src/Performance.cuh +++ b/engine/src/Performance.cuh @@ -32,6 +32,8 @@ namespace Performance { void test_with_other_solutions(const std::string& integral_str, PrintRoutine print_results = print_human_readable_results); + + void test_performance(const std::vector& integrals, PrintRoutine print_results = print_human_readable_results); } #endif \ No newline at end of file diff --git a/engine/src/main.cu b/engine/src/main.cu index e7dc672..6850847 100644 --- a/engine/src/main.cu +++ b/engine/src/main.cu @@ -43,17 +43,6 @@ std::string e_tower(size_t n) { return res; } -void warmup() { - Performance::test_with_other_solutions("x+sin(e^x)*e^x", Performance::do_not_print_results); -} - -void test_performance(const std::vector& integrals) { - for (auto integral_str : integrals) { - Performance::test_with_other_solutions(integral_str, - Performance::print_human_readable_results); - } -} - int main() { if constexpr (Consts::DEBUG) { fmt::print("Running in debug mode\n"); @@ -61,46 +50,48 @@ int main() { Sym::Static::init_functions(); - printf("Warmup...\n\n"); - warmup(); + fmt::print("First two integrals are for warmup!\n"); - test_performance({ - "x^4-5x+3", + Performance::test_performance({ + "x+sin(x)", + "e^x", + "x^4-5*x+3", "(1-x^3)^2", - "x^3t^2", + "x^3*t^2", "x^(2/3)", "1/x^(2/3)", - "(x^3-3x^2+1)/x^5", + "(x^3-3*x^2+1)/x^5", "(x^(2/3)-sqrt(x))/x^(1/3)", - "2e^x+3*5^x", - "x cos(x)", - "x e^x", - "x^2e^x", - "x^2cos(x)", - "sqrt(x)ln(x)", + "2*e^x+3*5^x", + "x*cos(x)", + "x*e^x", + "x^2*e^x", + "x^2*cos(x)", + "sqrt(x)*ln(x)", "ln(x)/x^4", - "(2x-1)^20", - "cos(3x-1)", + "(2*x-1)^20", + "cos(3*x-1)", "e^x/(3+e^x)", - "tg(x)", - "tg(x)/cos^2(x)", - "cos(x)e^sin(x)", + "tan(x)", + "tan(x)/cos(x)^2", + "cos(x)*e^sin(x)", "cos(x)/sqrt(1+sin(x))", "6^(1-x)", "(sqrt(x)-2)^2/x^2", - "1/(sin^2(x)cos^2(x))", - "x*arctg(x)", + "1/(sin(x)^2*cos(x)^2)", + "x*arctan(x)", "abs(x)", "(abs(1-x^2)+1+x^2)/2", - "tg^2(x)", - "(2x-3)^10", - "sin^5(x)cos(x)", + "tan(x)^2", + "(2*x-3)^10", + "sin(x)^5*cos(x)", "ln(x)", "x*cos(x)", - "x^2e^(1-x)", - "sin^5(x)", - "sin^4(x)cos^3(x)", + "x^2*e^(1-x)", + "sin(x)^5", + "sin(x)^4*cos(x)^3", "e^x*e^e^x*e^e^e^x*e^e^e^e^x*e^e^e^e^e^x*e^e^e^e^e^e^x", + "e^x*(x+1)*ln(x)", }); // const auto integral = Sym::integral(Parser::parse_function("tg(x)^2")); From 89cb3d00672d62f2f28a2fd045e9118e8e703d6c Mon Sep 17 00:00:00 2001 From: zh-qs Date: Thu, 19 Jan 2023 17:04:41 +0100 Subject: [PATCH 06/13] change timeit to tic-toc --- engine/src/Performance.cu | 8 +++++--- engine/src/main.cu | 17 ++++++++++++----- 2 files changed, 17 insertions(+), 8 deletions(-) diff --git a/engine/src/Performance.cu b/engine/src/Performance.cu index 064ee05..bdf426d 100644 --- a/engine/src/Performance.cu +++ b/engine/src/Performance.cu @@ -127,7 +127,8 @@ namespace Performance { replace(expression, "sgn", "sign"); replace(expression, "ln", "log"); - result += fmt::format(R"(f=@()int({},x);fprintf("%f\n",timeit(f));)", expression); + //result += fmt::format(R"(f=@()int({},x);fprintf("%f\n",timeit(f));)", expression); + result += fmt::format(R"(tic;int({},x);t=toc;fprintf("%f\n",t);)", expression); } return result + "\'"; @@ -190,9 +191,10 @@ namespace Performance { for (int i = 0; i < integrals.size(); ++i) { Sym::Integrator integrator; + const auto integral = Sym::integral(Parser::parse_function(integrals[i])); const clock_t start = clock(); const auto result = - integrator.solve_integral(Sym::integral(Parser::parse_function(integrals[i]))); + integrator.solve_integral(integral); const clock_t end = clock(); cusymint_seconds_vector[i] = static_cast(end - start) / CLOCKS_PER_SEC; @@ -202,7 +204,7 @@ namespace Performance { printf("Computing on Mathematica...\n"); auto mathematica_result = exec_and_read_output(make_mathematica_command_batch(integrals)); printf("Computing on SymPy...\n"); - auto sympy_result = exec_and_read_output(make_sympy_command_batch(integrals)); + auto sympy_result = std::string("");//exec_and_read_output(make_sympy_command_batch(integrals)); printf("Computing on MATLAB...\n"); auto matlab_result = exec_and_read_output(make_matlab_command_batch(integrals)); diff --git a/engine/src/main.cu b/engine/src/main.cu index 6850847..9117497 100644 --- a/engine/src/main.cu +++ b/engine/src/main.cu @@ -1,7 +1,9 @@ #include +#include #include #include +#include #include #include @@ -52,7 +54,7 @@ int main() { fmt::print("First two integrals are for warmup!\n"); - Performance::test_performance({ + std::vector integrals = { "x+sin(x)", "e^x", "x^4-5*x+3", @@ -92,16 +94,21 @@ int main() { "sin(x)^4*cos(x)^3", "e^x*e^e^x*e^e^e^x*e^e^e^e^x*e^e^e^e^e^x*e^e^e^e^e^e^x", "e^x*(x+1)*ln(x)", - }); + }; - // const auto integral = Sym::integral(Parser::parse_function("tg(x)^2")); + // add huge sum + // integrals.push_back(fmt::format("{}", fmt::join(integrals, "+"))); + + Performance::test_performance(integrals); + + // const auto integral = Sym::integral(Parser::parse_function(integrals.back())); // fmt::print("Trying to solve an integral: {}\n", integral.data()->to_tex()); // Sym::Integrator integrator; // Sym::ComputationHistory history; - // const auto solution = integrator.solve_integral_with_history(integral, history); - + // const auto solution = integrator.solve_integral(integral);//integrator.solve_integral_with_history(integral, history); + // if (solution.has_value()) { // fmt::print("Success! Solution:\n{} + C\n", solution.value().data()->to_tex()); From d41fcb9d0aa1287d16b285cf47e8051f7c17ec42 Mon Sep 17 00:00:00 2001 From: zh-qs Date: Fri, 20 Jan 2023 08:39:52 +0100 Subject: [PATCH 07/13] Add performance tests generating --- engine/src/Performance.cu | 34 +++++++++++++-- engine/src/Performance.cuh | 14 ++++-- engine/src/main.cu | 87 ++++++++++++++++++++++++++++++++++++-- 3 files changed, 124 insertions(+), 11 deletions(-) diff --git a/engine/src/Performance.cu b/engine/src/Performance.cu index bdf426d..db2af85 100644 --- a/engine/src/Performance.cu +++ b/engine/src/Performance.cu @@ -127,7 +127,8 @@ namespace Performance { replace(expression, "sgn", "sign"); replace(expression, "ln", "log"); - //result += fmt::format(R"(f=@()int({},x);fprintf("%f\n",timeit(f));)", expression); + // result += fmt::format(R"(f=@()int({},x);fprintf("%f\n",timeit(f));)", + // expression); result += fmt::format(R"(tic;int({},x);t=toc;fprintf("%f\n",t);)", expression); } @@ -161,6 +162,10 @@ namespace Performance { matlab_result.c_str()); } + void print_csv_headers() { + printf("integral;solved_by_cusymint;cusymint_time[s];mathematica_time[s];sympy_time[s];matlab_time[s]\n"); + } + void test_with_other_solutions(const std::string& integral_str, PrintRoutine print_results) { const auto integral = Sym::integral(Parser::parse_function(integral_str)); @@ -180,6 +185,24 @@ namespace Performance { sympy_result, matlab_result); } + void test_cuda_and_print_commands(const std::vector& integrals) { + printf("Computing on CUDA...\n"); + + for (int i = 0; i < integrals.size(); ++i) { + Sym::Integrator integrator; + const auto integral = Sym::integral(Parser::parse_function(integrals[i])); + const clock_t start = clock(); + const auto result = integrator.solve_integral(integral); + const clock_t end = clock(); + + fmt::print("{}{}\n", static_cast(end - start) / CLOCKS_PER_SEC, result.has_value() ? "" : " (failure)"); + } + + printf("%s\n", make_mathematica_command_batch(integrals).c_str()); + printf("%s\n", make_sympy_command_batch(integrals).c_str()); + printf("%s\n", make_matlab_command_batch(integrals).c_str()); + } + void test_performance(const std::vector& integrals, PrintRoutine print_results) { std::vector cusymint_seconds_vector(integrals.size()); std::vector cusymint_results_vector(integrals.size()); @@ -193,8 +216,7 @@ namespace Performance { Sym::Integrator integrator; const auto integral = Sym::integral(Parser::parse_function(integrals[i])); const clock_t start = clock(); - const auto result = - integrator.solve_integral(integral); + const auto result = integrator.solve_integral(integral); const clock_t end = clock(); cusymint_seconds_vector[i] = static_cast(end - start) / CLOCKS_PER_SEC; @@ -204,7 +226,7 @@ namespace Performance { printf("Computing on Mathematica...\n"); auto mathematica_result = exec_and_read_output(make_mathematica_command_batch(integrals)); printf("Computing on SymPy...\n"); - auto sympy_result = std::string("");//exec_and_read_output(make_sympy_command_batch(integrals)); + auto sympy_result = exec_and_read_output(make_sympy_command_batch(integrals)); printf("Computing on MATLAB...\n"); auto matlab_result = exec_and_read_output(make_matlab_command_batch(integrals)); @@ -231,6 +253,10 @@ namespace Performance { } } + if (print_results == print_csv_results) { + print_csv_headers(); + } + for (int i = 0; i < integrals.size(); ++i) { print_results(integrals[i], cusymint_seconds_vector[i], cusymint_results_vector[i], mathematica_results_vector[i], sympy_results_vector[i], diff --git a/engine/src/Performance.cuh b/engine/src/Performance.cuh index d22fcba..a609f4a 100644 --- a/engine/src/Performance.cuh +++ b/engine/src/Performance.cuh @@ -6,6 +6,7 @@ #include #include #include +#include #include "Evaluation/Integrator.cuh" #include "Parser/Parser.cuh" @@ -16,7 +17,8 @@ namespace Performance { using PrintRoutine = void (*)(const std::string& integral_str, const double& cusymint_seconds, bool cusymint_success, const std::string& mathematica_result, - const std::string& sympy_result, const std::string& matlab_result); + const std::string& sympy_result, + const std::string& matlab_result); void do_not_print_results(const std::string& integral_str, const double& cusymint_seconds, bool cusymint_success, const std::string& mathematica_result, @@ -24,16 +26,22 @@ namespace Performance { void print_human_readable_results(const std::string& integral_str, const double& cusymint_seconds, bool cusymint_success, const std::string& mathematica_result, - const std::string& sympy_result, const std::string& matlab_result); + const std::string& sympy_result, + const std::string& matlab_result); void print_csv_results(const std::string& integral_str, const double& cusymint_seconds, bool cusymint_success, const std::string& mathematica_result, const std::string& sympy_result, const std::string& matlab_result); + void print_csv_headers(); + void test_with_other_solutions(const std::string& integral_str, PrintRoutine print_results = print_human_readable_results); - void test_performance(const std::vector& integrals, PrintRoutine print_results = print_human_readable_results); + void test_cuda_and_print_commands(const std::vector& integrals); + + void test_performance(const std::vector& integrals, + PrintRoutine print_results = print_human_readable_results); } #endif \ No newline at end of file diff --git a/engine/src/main.cu b/engine/src/main.cu index 9117497..931636d 100644 --- a/engine/src/main.cu +++ b/engine/src/main.cu @@ -1,3 +1,4 @@ +#include #include #include @@ -45,6 +46,64 @@ std::string e_tower(size_t n) { return res; } +std::vector generate_e_towers(size_t end) { + std::vector result(end + 1); + result[0] = "1"; + if (end > 0) { + result[1] = "e^x"; + } + for (size_t i = 2; i <= end; ++i) { + std::string e_power = "*"; + for (size_t j = 1; j <= i; ++j) { + e_power += "e^"; + } + e_power += "x"; + result[i] = result[i - 1] + e_power; + } + return result; +} + +std::vector generate_sin_towers(size_t end) { + std::vector result(end + 1); + result[0] = "cos(x)"; + for (size_t i = 1; i <= end; ++i) { + std::string trig = "*cos("; + for (size_t j = 0; j < i; ++j) { + trig += "sin("; + } + trig += "x"; + for (size_t j = 0; j <= i; ++j) { + trig += ")"; + } + result[i] = result[i - 1] + trig; + } + return result; +} + +std::vector generate_geometric_sums(size_t start, size_t end, size_t step = 1) { + std::vector result; + for (size_t i = start; i <= end; i += step) { + result.push_back(fmt::format("(x^{}-1)/(x-1)", i)); + } + return result; +} + +std::vector generate_random_polynomials(size_t min_rank, size_t max_rank, size_t count, + size_t seed) { + constexpr size_t MAX_COEFF = 50; + std::vector result(count); + srand(seed); + for (size_t i = 0; i < count; ++i) { + size_t rank = (rand() % (max_rank - min_rank + 1)) + min_rank; + result[i] = fmt::format("{}", rand() % (MAX_COEFF + 1)); + for (size_t j = 1; j <= rank; ++j) { + result[i] += fmt::format("{}{}*x^{}", (rand() % 2 == 0) ? "+" : "-", + rand() % (MAX_COEFF + 1), j); + } + } + return result; +} + int main() { if constexpr (Consts::DEBUG) { fmt::print("Running in debug mode\n"); @@ -92,14 +151,33 @@ int main() { "x^2*e^(1-x)", "sin(x)^5", "sin(x)^4*cos(x)^3", - "e^x*e^e^x*e^e^e^x*e^e^e^e^x*e^e^e^e^e^x*e^e^e^e^e^e^x", "e^x*(x+1)*ln(x)", + "(1-x^100)/(1-x)", + "(1-sin(x)^10)/(1-sin(x))*cos(x)", }; // add huge sum // integrals.push_back(fmt::format("{}", fmt::join(integrals, "+"))); - Performance::test_performance(integrals); + //fmt::print("{}", fmt::join(generate_sin_towers(5), ",")); + + // add e towers + // const auto e_towers = generate_e_towers(10); + // integrals.insert(integrals.end(), e_towers.begin(), e_towers.end()); + + // add sin towers + // const auto sin_towers = generate_sin_towers(8); + // integrals.insert(integrals.end(), sin_towers.begin(), sin_towers.end()); + + // add geometric sums + //const auto geo_sums = generate_geometric_sums(0, 200, 5); + //integrals.insert(integrals.end(), geo_sums.begin(), geo_sums.end()); + + // add random polynomials + //const auto polynomials = generate_random_polynomials(100, 300, 10, 2137); + //integrals.insert(integrals.end(), polynomials.begin(), polynomials.end()); + + Performance::test_cuda_and_print_commands(integrals); // const auto integral = Sym::integral(Parser::parse_function(integrals.back())); @@ -107,8 +185,9 @@ int main() { // Sym::Integrator integrator; // Sym::ComputationHistory history; - // const auto solution = integrator.solve_integral(integral);//integrator.solve_integral_with_history(integral, history); - + // const auto solution = + // integrator.solve_integral_with_history(integral, history); + // if (solution.has_value()) { // fmt::print("Success! Solution:\n{} + C\n", solution.value().data()->to_tex()); From 62688ec02f207118e1995b845e185fc47b8ef796 Mon Sep 17 00:00:00 2001 From: zh-qs Date: Fri, 20 Jan 2023 09:48:37 +0100 Subject: [PATCH 08/13] Clear Mathematica cache --- engine/src/Performance.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/engine/src/Performance.cu b/engine/src/Performance.cu index db2af85..dd8d39c 100644 --- a/engine/src/Performance.cu +++ b/engine/src/Performance.cu @@ -91,7 +91,7 @@ namespace Performance { capitalize_substring(expression, "pi"); capitalize_substring(expression, "sqrt"); result += fmt::format("{{t,b}}=AbsoluteTiming[Integrate[ToExpression[\"{}\"," - "TraditionalForm],x]];Print[t];", + "TraditionalForm],x]];Print[t];ClearSystemCache[];", expression); } From 76f4222f41edef7f89fcf9df03db719bafd56852 Mon Sep 17 00:00:00 2001 From: zh-qs Date: Fri, 20 Jan 2023 10:16:37 +0100 Subject: [PATCH 09/13] add warning --- engine/src/Performance.cu | 1 + 1 file changed, 1 insertion(+) diff --git a/engine/src/Performance.cu b/engine/src/Performance.cu index dd8d39c..4956e94 100644 --- a/engine/src/Performance.cu +++ b/engine/src/Performance.cu @@ -225,6 +225,7 @@ namespace Performance { printf("Computing on Mathematica...\n"); auto mathematica_result = exec_and_read_output(make_mathematica_command_batch(integrals)); + // Warning: computing na integral which requires many substitutions on SymPy may hang your computer and is extremely slow! printf("Computing on SymPy...\n"); auto sympy_result = exec_and_read_output(make_sympy_command_batch(integrals)); printf("Computing on MATLAB...\n"); From 5d20c875fd6b1980d2b7aee5119456d047345287 Mon Sep 17 00:00:00 2001 From: zh-qs Date: Sat, 21 Jan 2023 01:29:17 +0100 Subject: [PATCH 10/13] Add memory testing --- engine/src/Evaluation/Integrator.cu | 5 +- engine/src/Evaluation/Integrator.cuh | 125 +++++++++++++++++++++++++++ engine/src/Performance.cu | 36 +++++++- engine/src/Performance.cuh | 2 + engine/src/main.cu | 111 ++++++++++++++---------- 5 files changed, 227 insertions(+), 52 deletions(-) diff --git a/engine/src/Evaluation/Integrator.cu b/engine/src/Evaluation/Integrator.cu index e33b3f1..5998392 100644 --- a/engine/src/Evaluation/Integrator.cu +++ b/engine/src/Evaluation/Integrator.cu @@ -4,9 +4,9 @@ #include #include +#include "ComputationHistory.cuh" #include "IntegratorKernels.cuh" #include "StaticFunctions.cuh" -#include "ComputationHistory.cuh" #include "Evaluation/Heuristic/Heuristic.cuh" #include "Symbol/SubexpressionCandidate.cuh" @@ -188,7 +188,8 @@ namespace Sym { scan_array_2.zero_mem(); cudaDeviceSynchronize(); - help_space.reoffset_like(integrals.iterator(), HELP_SPACE_MULTIPLIER, Heuristic::COUNT); + help_space.reoffset_like(integrals.iterator(), + HELP_SPACE_MULTIPLIER, Heuristic::COUNT); Kernel::check_heuristics_applicability<<>>( integrals, expressions, help_space, scan_array_1, scan_array_2); diff --git a/engine/src/Evaluation/Integrator.cuh b/engine/src/Evaluation/Integrator.cuh index ed3077a..a99c74d 100644 --- a/engine/src/Evaluation/Integrator.cuh +++ b/engine/src/Evaluation/Integrator.cuh @@ -238,6 +238,131 @@ namespace Sym { return std::nullopt; } + + size_t memory_usage_for_integral(const std::vector& integral, + size_t& total_memory) { + size_t initial_usage; + size_t max_usage; + size_t usage; + + size_t initial_size = + sizeof(Symbol) * + (expressions.symbols_capacity() + integrals.symbols_capacity() + + expressions_swap.symbols_capacity() + integrals_swap.symbols_capacity() + + help_space.symbols_capacity()) + + sizeof(size_t) * (scan_array_1.size() + scan_array_2.size()) + + sizeof(EvaluationStatus) * + (evaluation_statuses_1.size() + evaluation_statuses_2.size()); + size_t max_shown_size = initial_size; + + cudaMemGetInfo(&initial_usage, &total_memory); + max_usage = initial_usage; + + expressions.load_from_vector({single_integral_vacancy()}); + integrals.load_from_vector({first_expression_candidate(integral)}); + + max_shown_size = std::max( + max_shown_size, + sizeof(Symbol) * + (expressions.symbols_capacity() + integrals.symbols_capacity() + + expressions_swap.symbols_capacity() + integrals_swap.symbols_capacity() + + help_space.symbols_capacity()) + + sizeof(size_t) * (scan_array_1.size() + scan_array_2.size()) + + sizeof(EvaluationStatus) * + (evaluation_statuses_1.size() + evaluation_statuses_2.size())); + + cudaMemGetInfo(&usage, nullptr); + max_usage = std::min(usage, max_usage); + + for (size_t i = 0;; ++i) { + simplify_integrals(); + + cudaMemGetInfo(&usage, nullptr); + max_usage = std::min(usage, max_usage); + + max_shown_size = std::max( + max_shown_size, + sizeof(Symbol) * + (expressions.symbols_capacity() + integrals.symbols_capacity() + + expressions_swap.symbols_capacity() + + integrals_swap.symbols_capacity() + help_space.symbols_capacity()) + + sizeof(size_t) * (scan_array_1.size() + scan_array_2.size()) + + sizeof(EvaluationStatus) * + (evaluation_statuses_1.size() + evaluation_statuses_2.size())); + + check_for_known_integrals(); + apply_known_integrals(); + + cudaMemGetInfo(&usage, nullptr); + max_usage = std::min(usage, max_usage); + + max_shown_size = std::max( + max_shown_size, + sizeof(Symbol) * + (expressions.symbols_capacity() + integrals.symbols_capacity() + + expressions_swap.symbols_capacity() + + integrals_swap.symbols_capacity() + help_space.symbols_capacity()) + + sizeof(size_t) * (scan_array_1.size() + scan_array_2.size()) + + sizeof(EvaluationStatus) * + (evaluation_statuses_1.size() + evaluation_statuses_2.size())); + + if (is_original_expression_solved()) { + break; + } + + remove_unnecessary_candidates(); + + cudaMemGetInfo(&usage, nullptr); + max_usage = std::min(usage, max_usage); + + max_shown_size = std::max( + max_shown_size, + sizeof(Symbol) * + (expressions.symbols_capacity() + integrals.symbols_capacity() + + expressions_swap.symbols_capacity() + + integrals_swap.symbols_capacity() + help_space.symbols_capacity()) + + sizeof(size_t) * (scan_array_1.size() + scan_array_2.size()) + + sizeof(EvaluationStatus) * + (evaluation_statuses_1.size() + evaluation_statuses_2.size())); + + check_heuristics_applicability(); + apply_heuristics(); + + cudaMemGetInfo(&usage, nullptr); + max_usage = std::min(usage, max_usage); + + max_shown_size = std::max( + max_shown_size, + sizeof(Symbol) * + (expressions.symbols_capacity() + integrals.symbols_capacity() + + expressions_swap.symbols_capacity() + + integrals_swap.symbols_capacity() + help_space.symbols_capacity()) + + sizeof(size_t) * (scan_array_1.size() + scan_array_2.size()) + + sizeof(EvaluationStatus) * + (evaluation_statuses_1.size() + evaluation_statuses_2.size())); + + if (has_original_expression_failed()) { + break; + } + + remove_failed_candidates(); + + cudaMemGetInfo(&usage, nullptr); + max_usage = std::min(usage, max_usage); + + max_shown_size = std::max( + max_shown_size, + sizeof(Symbol) * + (expressions.symbols_capacity() + integrals.symbols_capacity() + + expressions_swap.symbols_capacity() + + integrals_swap.symbols_capacity() + help_space.symbols_capacity()) + + sizeof(size_t) * (scan_array_1.size() + scan_array_2.size()) + + sizeof(EvaluationStatus) * + (evaluation_statuses_1.size() + evaluation_statuses_2.size())); + } + printf("Declared bytes:%lu\n", max_shown_size); + return initial_usage - max_usage; + } }; } diff --git a/engine/src/Performance.cu b/engine/src/Performance.cu index 4956e94..a12447b 100644 --- a/engine/src/Performance.cu +++ b/engine/src/Performance.cu @@ -1,5 +1,7 @@ #include "Performance.cuh" #include "Symbol/Integral.cuh" +#include "Symbol/SymbolType.cuh" +#include #include #include #include @@ -163,7 +165,8 @@ namespace Performance { } void print_csv_headers() { - printf("integral;solved_by_cusymint;cusymint_time[s];mathematica_time[s];sympy_time[s];matlab_time[s]\n"); + printf("integral;solved_by_cusymint;cusymint_time[s];mathematica_time[s];sympy_time[s];" + "matlab_time[s]\n"); } void test_with_other_solutions(const std::string& integral_str, PrintRoutine print_results) { @@ -195,7 +198,8 @@ namespace Performance { const auto result = integrator.solve_integral(integral); const clock_t end = clock(); - fmt::print("{}{}\n", static_cast(end - start) / CLOCKS_PER_SEC, result.has_value() ? "" : " (failure)"); + fmt::print("{}{}\n", static_cast(end - start) / CLOCKS_PER_SEC, + result.has_value() ? "" : " (failure)"); } printf("%s\n", make_mathematica_command_batch(integrals).c_str()); @@ -203,6 +207,31 @@ namespace Performance { printf("%s\n", make_matlab_command_batch(integrals).c_str()); } + void test_memory_occupance(const std::vector& integrals) { + printf("Memory occupance:\n"); + for (int i = 0; i < integrals.size(); ++i) { + printf(" Integral %s:\n", integrals[i].c_str()); + + size_t initial_mem; + size_t after_init; + + cudaMemGetInfo(&initial_mem, nullptr); + + Sym::Integrator integrator; + + cudaMemGetInfo(&after_init, nullptr); + initial_mem -= after_init; + + size_t total; + const auto integral = Sym::integral(Parser::parse_function(integrals[i])); + + const size_t usage = integrator.memory_usage_for_integral(integral, total); + + printf(" %10lu/%10lu B used (%f%%)\n", usage + initial_mem, total, + static_cast(usage + initial_mem) / total * 100); + } + } + void test_performance(const std::vector& integrals, PrintRoutine print_results) { std::vector cusymint_seconds_vector(integrals.size()); std::vector cusymint_results_vector(integrals.size()); @@ -225,7 +254,8 @@ namespace Performance { printf("Computing on Mathematica...\n"); auto mathematica_result = exec_and_read_output(make_mathematica_command_batch(integrals)); - // Warning: computing na integral which requires many substitutions on SymPy may hang your computer and is extremely slow! + // Warning: computing na integral which requires many substitutions on SymPy may hang your + // computer and is extremely slow! printf("Computing on SymPy...\n"); auto sympy_result = exec_and_read_output(make_sympy_command_batch(integrals)); printf("Computing on MATLAB...\n"); diff --git a/engine/src/Performance.cuh b/engine/src/Performance.cuh index a609f4a..ced99b5 100644 --- a/engine/src/Performance.cuh +++ b/engine/src/Performance.cuh @@ -38,6 +38,8 @@ namespace Performance { void test_with_other_solutions(const std::string& integral_str, PrintRoutine print_results = print_human_readable_results); + void test_memory_occupance(const std::vector& integrals); + void test_cuda_and_print_commands(const std::vector& integrals); void test_performance(const std::vector& integrals, diff --git a/engine/src/main.cu b/engine/src/main.cu index 931636d..369883b 100644 --- a/engine/src/main.cu +++ b/engine/src/main.cu @@ -104,6 +104,22 @@ std::vector generate_random_polynomials(size_t min_rank, size_t max return result; } +std::vector generate_random_trig_polynomials(size_t start, size_t end, size_t step, + size_t seed) { + constexpr size_t MAX_COEFF = 50; + std::vector result; + srand(seed); + for (size_t rank = start; rank <= end; rank += step) { + std::string str = fmt::format("{}", rand() % (MAX_COEFF + 1)); + for (size_t j = 1; j <= rank; ++j) { + str += fmt::format("{}{}*sin(x)^{}*cos(x)", (rand() % 2 == 0) ? "+" : "-", + rand() % (MAX_COEFF + 1), j); + } + result.push_back(str); + } + return result; +} + int main() { if constexpr (Consts::DEBUG) { fmt::print("Running in debug mode\n"); @@ -114,52 +130,53 @@ int main() { fmt::print("First two integrals are for warmup!\n"); std::vector integrals = { - "x+sin(x)", - "e^x", - "x^4-5*x+3", - "(1-x^3)^2", - "x^3*t^2", - "x^(2/3)", - "1/x^(2/3)", - "(x^3-3*x^2+1)/x^5", - "(x^(2/3)-sqrt(x))/x^(1/3)", - "2*e^x+3*5^x", - "x*cos(x)", - "x*e^x", - "x^2*e^x", - "x^2*cos(x)", - "sqrt(x)*ln(x)", - "ln(x)/x^4", - "(2*x-1)^20", - "cos(3*x-1)", - "e^x/(3+e^x)", - "tan(x)", - "tan(x)/cos(x)^2", - "cos(x)*e^sin(x)", - "cos(x)/sqrt(1+sin(x))", - "6^(1-x)", - "(sqrt(x)-2)^2/x^2", - "1/(sin(x)^2*cos(x)^2)", - "x*arctan(x)", - "abs(x)", - "(abs(1-x^2)+1+x^2)/2", - "tan(x)^2", - "(2*x-3)^10", - "sin(x)^5*cos(x)", - "ln(x)", - "x*cos(x)", - "x^2*e^(1-x)", - "sin(x)^5", - "sin(x)^4*cos(x)^3", - "e^x*(x+1)*ln(x)", - "(1-x^100)/(1-x)", - "(1-sin(x)^10)/(1-sin(x))*cos(x)", + "x+sin(x)", "e^x", + // "x^4-5*x+3", + // "(1-x^3)^2", + // "x^3*t^2", + // "x^(2/3)", + // "1/x^(2/3)", + // "(x^3-3*x^2+1)/x^5", + // "(x^(2/3)-sqrt(x))/x^(1/3)", + // "2*e^x+3*5^x", + // "x*cos(x)", + // "x*e^x", + // "x^2*e^x", + // "x*e^(-x)", + // "x^2*cos(x)", + // "sqrt(x)*ln(x)", + // "ln(x)/x^4", + // "(2*x-1)^20", + // "cos(3*x-1)", + // "e^x/(3+e^x)", + // "tan(x)", + // "tan(x)/cos(x)^2", + // "cos(x)*e^sin(x)", + // "cos(x)/sqrt(1+sin(x))", + // "6^(1-x)", + // "(sqrt(x)-2)^2/x^2", + // "1/(sin(x)^2*cos(x)^2)", + // "x*arctan(x)", + // "abs(x)", + // "(abs(1-x^2)+1+x^2)/2", + // "tan(x)^2", + // "(2*x-3)^10", + // "sin(x)^5*cos(x)", + // "ln(x)", + // "x*cos(x)", + // "x^2*e^(1-x)", + // "sin(x)^5", + // "sin(x)^4*cos(x)^3", + // "e^x*(x+1)*ln(x)", + // "(1-x^200)/(1-x)", + // "(1-sin(x)^20)/(1-sin(x))*cos(x)", + // e_tower(70), }; // add huge sum // integrals.push_back(fmt::format("{}", fmt::join(integrals, "+"))); - //fmt::print("{}", fmt::join(generate_sin_towers(5), ",")); + // fmt::print("{}", fmt::join(generate_sin_towers(5), ",")); // add e towers // const auto e_towers = generate_e_towers(10); @@ -170,15 +187,15 @@ int main() { // integrals.insert(integrals.end(), sin_towers.begin(), sin_towers.end()); // add geometric sums - //const auto geo_sums = generate_geometric_sums(0, 200, 5); - //integrals.insert(integrals.end(), geo_sums.begin(), geo_sums.end()); + const auto geo_sums = generate_geometric_sums(0, 200, 5); + integrals.insert(integrals.end(), geo_sums.begin(), geo_sums.end()); // add random polynomials - //const auto polynomials = generate_random_polynomials(100, 300, 10, 2137); - //integrals.insert(integrals.end(), polynomials.begin(), polynomials.end()); - - Performance::test_cuda_and_print_commands(integrals); + const auto polynomials = generate_random_trig_polynomials(0, 5, 100, 2137); + integrals.insert(integrals.end(), polynomials.begin(), polynomials.end()); + // Performance::test_performance(integrals, Performance::print_csv_results); + Performance::test_memory_occupance(integrals); // const auto integral = Sym::integral(Parser::parse_function(integrals.back())); // fmt::print("Trying to solve an integral: {}\n", integral.data()->to_tex()); From 5b9338a71fbf512d82f4d9a4c3de4821fd45cc92 Mon Sep 17 00:00:00 2001 From: zh-qs Date: Sat, 21 Jan 2023 21:15:50 +0100 Subject: [PATCH 11/13] Add memory occupance benchmark for history --- engine/src/Evaluation/ComputationHistory.cu | 16 ++++++++++++++++ engine/src/Evaluation/ComputationHistory.cuh | 3 +++ engine/src/Performance.cu | 15 +++++++++++++++ engine/src/Performance.cuh | 2 ++ engine/src/main.cu | 5 +++-- 5 files changed, 39 insertions(+), 2 deletions(-) diff --git a/engine/src/Evaluation/ComputationHistory.cu b/engine/src/Evaluation/ComputationHistory.cu index 9118ed4..242c5fe 100644 --- a/engine/src/Evaluation/ComputationHistory.cu +++ b/engine/src/Evaluation/ComputationHistory.cu @@ -234,6 +234,22 @@ namespace Sym { fmt::print("}}\n"); } + size_t ComputationStep::get_memory_occupance_in_bytes() const { + size_t size = 0; + for (const auto& expr : expression_tree) { + size += expr.size(); + } + return sizeof(ComputationStep) + size * sizeof(Symbol); + } + + size_t ComputationHistory::get_memory_occupance_in_bytes() const { + size_t size = 0; + for (const auto& step : computation_steps) { + size += step.get_memory_occupance_in_bytes(); + } + return sizeof(ComputationHistory) + size; + } + void ComputationHistory::complete() { if (!is_solution_at_end_of_steps()) { Util::crash("ComputationHistory must have solution at the end to complete"); diff --git a/engine/src/Evaluation/ComputationHistory.cuh b/engine/src/Evaluation/ComputationHistory.cuh index 01b283a..474c533 100644 --- a/engine/src/Evaluation/ComputationHistory.cuh +++ b/engine/src/Evaluation/ComputationHistory.cuh @@ -33,6 +33,7 @@ namespace Sym { std::vector get_expression() const; TransformationList get_operations(const ComputationStep& previous_step) const; void print_step() const; + size_t get_memory_occupance_in_bytes() const; }; using ComputationStepCollection = std::list; @@ -58,6 +59,8 @@ namespace Sym { void print_history() const; const ComputationStepCollection& get_steps() const { return computation_steps; } + + size_t get_memory_occupance_in_bytes() const; }; } diff --git a/engine/src/Performance.cu b/engine/src/Performance.cu index a12447b..2d356fb 100644 --- a/engine/src/Performance.cu +++ b/engine/src/Performance.cu @@ -232,6 +232,21 @@ namespace Performance { } } + void test_history_cpu_memory_occupance(const std::vector& integrals) { + printf("Memory occupance for history:\n"); + for (int i = 0; i < integrals.size(); ++i) { + printf(" Integral %s:\n", integrals[i].c_str()); + + Sym::Integrator integrator; + Sym::ComputationHistory history; + const auto integral = Sym::integral(Parser::parse_function(integrals[i])); + + const auto result = integrator.solve_integral_with_history(integral, history); + + printf(" %10luB used%s\n", history.get_memory_occupance_in_bytes(), result.has_value() ? "" : " (failure)"); + } + } + void test_performance(const std::vector& integrals, PrintRoutine print_results) { std::vector cusymint_seconds_vector(integrals.size()); std::vector cusymint_results_vector(integrals.size()); diff --git a/engine/src/Performance.cuh b/engine/src/Performance.cuh index ced99b5..85f5a84 100644 --- a/engine/src/Performance.cuh +++ b/engine/src/Performance.cuh @@ -40,6 +40,8 @@ namespace Performance { void test_memory_occupance(const std::vector& integrals); + void test_history_cpu_memory_occupance(const std::vector& integrals); + void test_cuda_and_print_commands(const std::vector& integrals); void test_performance(const std::vector& integrals, diff --git a/engine/src/main.cu b/engine/src/main.cu index 369883b..103403a 100644 --- a/engine/src/main.cu +++ b/engine/src/main.cu @@ -192,10 +192,11 @@ int main() { // add random polynomials const auto polynomials = generate_random_trig_polynomials(0, 5, 100, 2137); - integrals.insert(integrals.end(), polynomials.begin(), polynomials.end()); + //integrals.insert(integrals.end(), polynomials.begin(), polynomials.end()); // Performance::test_performance(integrals, Performance::print_csv_results); - Performance::test_memory_occupance(integrals); + // Performance::test_memory_occupance(integrals); + Performance::test_history_cpu_memory_occupance(integrals); // const auto integral = Sym::integral(Parser::parse_function(integrals.back())); // fmt::print("Trying to solve an integral: {}\n", integral.data()->to_tex()); From 5379dfe7a4424061ee1dc7be37268c2be46b04fd Mon Sep 17 00:00:00 2001 From: zh-qs Date: Sun, 22 Jan 2023 15:59:30 +0100 Subject: [PATCH 12/13] Make polynomial sequence not-random --- engine/src/main.cu | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/engine/src/main.cu b/engine/src/main.cu index 103403a..3b56373 100644 --- a/engine/src/main.cu +++ b/engine/src/main.cu @@ -104,16 +104,12 @@ std::vector generate_random_polynomials(size_t min_rank, size_t max return result; } -std::vector generate_random_trig_polynomials(size_t start, size_t end, size_t step, - size_t seed) { - constexpr size_t MAX_COEFF = 50; +std::vector generate_trig_polynomials(size_t start, size_t end, size_t step) { std::vector result; - srand(seed); for (size_t rank = start; rank <= end; rank += step) { - std::string str = fmt::format("{}", rand() % (MAX_COEFF + 1)); + std::string str = "1"; for (size_t j = 1; j <= rank; ++j) { - str += fmt::format("{}{}*sin(x)^{}*cos(x)", (rand() % 2 == 0) ? "+" : "-", - rand() % (MAX_COEFF + 1), j); + str += fmt::format("+sin(x)^{}*cos(x)", j); } result.push_back(str); } @@ -188,15 +184,15 @@ int main() { // add geometric sums const auto geo_sums = generate_geometric_sums(0, 200, 5); - integrals.insert(integrals.end(), geo_sums.begin(), geo_sums.end()); + //integrals.insert(integrals.end(), geo_sums.begin(), geo_sums.end()); // add random polynomials - const auto polynomials = generate_random_trig_polynomials(0, 5, 100, 2137); - //integrals.insert(integrals.end(), polynomials.begin(), polynomials.end()); + const auto polynomials = generate_trig_polynomials(0, 100, 5); + integrals.insert(integrals.end(), polynomials.begin(), polynomials.end()); // Performance::test_performance(integrals, Performance::print_csv_results); - // Performance::test_memory_occupance(integrals); - Performance::test_history_cpu_memory_occupance(integrals); + Performance::test_memory_occupance(integrals); + // Performance::test_history_cpu_memory_occupance(integrals); // const auto integral = Sym::integral(Parser::parse_function(integrals.back())); // fmt::print("Trying to solve an integral: {}\n", integral.data()->to_tex()); From ee0943a00f9e465ad8838456078c3c2dc3933ce4 Mon Sep 17 00:00:00 2001 From: zh-qs Date: Sun, 22 Jan 2023 16:12:09 +0100 Subject: [PATCH 13/13] Change 1 to cosine --- engine/src/main.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/engine/src/main.cu b/engine/src/main.cu index 3b56373..c5aedb2 100644 --- a/engine/src/main.cu +++ b/engine/src/main.cu @@ -107,7 +107,7 @@ std::vector generate_random_polynomials(size_t min_rank, size_t max std::vector generate_trig_polynomials(size_t start, size_t end, size_t step) { std::vector result; for (size_t rank = start; rank <= end; rank += step) { - std::string str = "1"; + std::string str = "cos(x)"; for (size_t j = 1; j <= rank; ++j) { str += fmt::format("+sin(x)^{}*cos(x)", j); }