diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index 983ffd8..8f0534e 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -1,4 +1,4 @@ -name: Python package +name: Python and C++ package on: pull_request: @@ -24,23 +24,39 @@ jobs: uses: actions/checkout@v2 with: submodules: recursive + - name: Set up Python ${{ matrix.python-version }} uses: actions/setup-python@v2 with: python-version: ${{ matrix.python-version }} + - name: Install dependencies run: | - python -m pip install --upgrade pip - pip install -r requirements.txt - pip install . + python3 -m pip install --upgrade pip + pip3 install -r requirements.txt + - name: Install testing packages run: | - # pip install flake8 - # pip install pydocstyle - pip install pytest - pip install pytest-cov pytest-xdist codecov - - name: Test with pytest + pip3 install pytest + pip3 install pytest-cov pytest-xdist codecov + + + # Below are the added steps for CMake + + - name: Setup CMake + uses: jwlawson/actions-setup-cmake@v1.9 + with: + cmake-version: 3.21.1 + + - name: Configure and run CMake + shell: bash run: | - pytest --cov=pyoscode tests + mkdir test_build + cd test_build + PYTHON_PATH=$(which python) + cmake .. -DCMAKE_BUILD_TYPE=Debug -DPYTHON_EXECUTABLE=${PYTHON_PATH} + cmake --build . --target test + - name: Upload coverage to codecov uses: codecov/codecov-action@v2 + diff --git a/.gitignore b/.gitignore index 243aece..0c922f2 100644 --- a/.gitignore +++ b/.gitignore @@ -4,10 +4,13 @@ !/JOSS-paper/* /examples/* !/examples/*.cpp +/examples/*.txt +examples/output.txt .eggs/* .vscode/* bin/* build/* +test_build/* dist/* env/* lib/* diff --git a/.gitmodules b/.gitmodules deleted file mode 100644 index 89d2db6..0000000 --- a/.gitmodules +++ /dev/null @@ -1,4 +0,0 @@ -[submodule "deps/eigen"] - path = deps/eigen - url = https://gitlab.com/libeigen/eigen - branch = 3.4 diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..cdedb33 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,140 @@ +cmake_minimum_required(VERSION 3.20.2) +project( + oscode + VERSION 1.3.0 + LANGUAGES C CXX) + +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED YES) +if(MSVC) + add_compile_options(/std:c++17) +endif() +set(CMAKE_CXX_EXTENSIONS NO) +if (CMAKE_BUILD_TYPE MATCHES Release) + set(CMAKE_VERBOSE_MAKEFILE NO) +else() + set(CMAKE_VERBOSE_MAKEFILE YES) +endif() + +# Build Types +set(CMAKE_BUILD_TYPE ${CMAKE_BUILD_TYPE} + CACHE STRING "Choose the type of build, options are: None Debug Release RelWithDebInfo MinSizeRel tsan asan lsan msan ubsan" + FORCE) + +# ThreadSanitizer +set(CMAKE_C_FLAGS_TSAN + "-fsanitize=thread -g -O1 -march=native -mtune=native" + CACHE STRING "Flags used by the C compiler during ThreadSanitizer builds." + FORCE) +set(CMAKE_CXX_FLAGS_TSAN + "-fsanitize=thread -g -O1 -march=native -mtune=native" + CACHE STRING "Flags used by the C++ compiler during ThreadSanitizer builds." + FORCE) + +# AddressSanitize +set(CMAKE_C_FLAGS_ASAN + "-fsanitize=address -fno-optimize-sibling-calls -fsanitize-address-use-after-scope -fno-omit-frame-pointer -g -Og -march=native -mtune=native" + CACHE STRING "Flags used by the C compiler during AddressSanitizer builds." + FORCE) +set(CMAKE_CXX_FLAGS_ASAN + "-fsanitize=address -fno-optimize-sibling-calls -fsanitize-address-use-after-scope -fno-omit-frame-pointer -g -Og -Wall -march=native -mtune=native" + CACHE STRING "Flags used by the C++ compiler during AddressSanitizer builds." + FORCE) + +# LeakSanitizer +set(CMAKE_C_FLAGS_LSAN + "-fsanitize=leak -fno-omit-frame-pointer -g -O1" + CACHE STRING "Flags used by the C compiler during LeakSanitizer builds." + FORCE) +set(CMAKE_CXX_FLAGS_LSAN + "-fsanitize=leak -fno-omit-frame-pointer -g -O1" + CACHE STRING "Flags used by the C++ compiler during LeakSanitizer builds." + FORCE) + + + +set(CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH}) +include(FetchContent) +# Externally provided libraries +FetchContent_Declare(googletest + GIT_REPOSITORY https://github.com/google/googletest.git + GIT_TAG main) + +FetchContent_MakeAvailable(googletest) + +FetchContent_Declare( + Eigen + GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git + GIT_TAG 3.4.0) +set(EIGEN_BUILD_DOC OFF) +# note: To disable eigen tests, +# you should put this code in a add_subdirectory to avoid to change +# BUILD_TESTING for your own project too since variables are directory +# scoped +set(BUILD_TESTING OFF) +set(EIGEN_BUILD_PKGCONFIG OFF) +FetchContent_MakeAvailable(Eigen) + +set(Boost_USE_STATIC_LIBS ON) +set(Boost_USE_STATIC_RUNTIME ON) +set(BOOST_ENABLE_CMAKE ON) +FetchContent_Declare( + boost + URL https://boostorg.jfrog.io/artifactory/main/release/1.82.0/source/boost_1_82_0.tar.gz +) + +FetchContent_GetProperties(boost) +if(NOT boost_POPULATED) + FetchContent_Populate(boost) + set(boost_INCLUDE_DIR ${boost_SOURCE_DIR}) +endif() + + +if (NOT CMAKE_CXX_COMPILER_ID STREQUAL "MSVC") + add_compile_options( + -DNO_FPRINTF_OUTPUT + -Wall + -O3) +endif() + +set(oscode_INCLUDE_DIR "${CMAKE_CURRENT_SOURCE_DIR}/include") + +set(PYTHON_EXECUTABLE "" CACHE FILEPATH "Python interpreter to use") + +option(OSCODE_HEADER_ONLY "Whether the headers are only included or if the python package should be built. Default builds the python package" OFF) + +if (NOT OSCODE_HEADER_ONLY) + if(NOT PYTHON_EXECUTABLE) + message(FATAL_ERROR "PYTHON_EXECUTABLE must be provided") + endif() + # Set the location of the eigen headers + #set(ENV{SOURCE_DIR} ${eigen_SOURCE_DIR}) + if (CMAKE_BUILD_TYPE MATCHES Release) + add_custom_target(pyoscode ALL + COMMAND ${CMAKE_COMMAND} -E env "OSCODE_EIGEN_INCLUDE_DIR=${eigen_SOURCE_DIR}" ${PYTHON_EXECUTABLE} setup.py build + COMMAND ${CMAKE_COMMAND} -E env "OSCODE_EIGEN_INCLUDE_DIR=${eigen_SOURCE_DIR}" pip3 install -e . + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} + COMMENT "Building Python package using setup.py" + ) + else() + add_custom_target(pyoscode ALL + COMMAND ${CMAKE_COMMAND} -E env "OSCODE_EIGEN_INCLUDE_DIR=${eigen_SOURCE_DIR}" ${PYTHON_EXECUTABLE} setup.py build_ext -i + COMMAND ${CMAKE_COMMAND} -E env "OSCODE_EIGEN_INCLUDE_DIR=${eigen_SOURCE_DIR}" pip3 install -e . + COMMAND ${CMAKE_COMMAND} -E env "OSCODE_EIGEN_INCLUDE_DIR=${eigen_SOURCE_DIR}" ${PYTHON_EXECUTABLE} setup.py build_ext -if + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} + COMMENT "Building Python package using setup.py" + ) + endif() +endif() +unset(OSCODE_HEADER_ONLY CACHE) +# Allows upstream users to call fetch_content etc. +add_library(oscode INTERFACE) + +target_include_directories(oscode PUBLIC INTERFACE + ${oscode_INCLUDE_DIR} +) + +target_link_libraries(oscode INTERFACE Eigen3::Eigen) + +include_directories(tests) +add_subdirectory(tests) diff --git a/LICENSE b/LICENSE index ca66dc7..456998e 100644 --- a/LICENSE +++ b/LICENSE @@ -1,6 +1,6 @@ BSD 3-Clause License -Copyright (c) 2020, Fruzsina Julia Agocs +Copyright (c) 2023, Fruzsina Julia Agocs All rights reserved. Redistribution and use in source and binary forms, with or without diff --git a/README.rst b/README.rst index cafa5b5..ed26b69 100644 --- a/README.rst +++ b/README.rst @@ -99,32 +99,27 @@ Python pip install pyoscode -or via the setup.py +Installing from source can be done via cmake with the following .. code:: bash - - git clone --recursive https://github.com/fruzsinaagocs/oscode - cd oscode - python setup.py install --user - -or - -.. code:: bash - - git clone --recursive https://github.com/fruzsinaagocs/oscode + git clone https://github.com/fruzsinaagocs/oscode cd oscode - pip install . + # For mac this will be $(which python3) + PYTHON_PATH=$(which python) + cmake -S . -B "build" -DCMAKE_BUILD_TYPE=Release -DPYTHON_EXECUTABLE=${PYTHON_PATH} + cd build + make oscode You can then import ``pyoscode`` from anywhere. Omit the ``--user`` option if you wish to install globally or in a virtual environment. If you have any difficulties, check out the `FAQs - Installation `__ section below. -You can check that things are working by running `tests/` (also ran by Travis continuous integration): +You can check that things are working by running the tests in `tests/` (also ran by github actions): .. code:: bash - pytest tests/ + make tests C++ ~~~ @@ -133,7 +128,7 @@ C++ .. code:: bash - git clone --recursive https://github.com/fruzsinaagocs/oscode + git clone https://github.com/fruzsinaagocs/oscode and then include the relevant header files in your C++ code: @@ -176,16 +171,18 @@ C++ :Introduction to oscode: `examples/burst.cpp` :To plot results from `burst.cpp`: `examples/plot_burst.py` -To compile and run: +To compile and run the burst example you can use the following: .. code:: bash cd examples/ - g++ -I../include/ -g -Wall -std=c++11 -c -o burst.o burst.cpp - g++ -I../include/ -g -Wall -std=c++11 -o burst burst.o - ./burst + cmake -S . -B "build" -DCMAKE_BUILD_TYPE=Release -DOSCODE_HEADER_ONLY=YES + cd build && make burst_ex && ./burst_ex +The `OSCODE_HEADER_ONLY` flag tells cmake to not create the python package target. +You can use the `CMakeLists.txt` file in the examples folder as a reference point for including `oscode` into your own package. + Documentation ------------- @@ -220,30 +217,41 @@ Further help You can get help by submitting an issue or posting a message on `Gitter `__. -FAQs ----- -Installation -~~~~~~~~~~~~ +Development +------------- -1. Eigen import errors: - .. code:: bash +Run the following to setup an environment for development: - pyoscode/_pyoscode.hpp:6:10: fatal error: Eigen/Dense: No such file or directory - #include - ^~~~~~~~~~~~~ + .. code:: bash + python -m venv env + source ./env/bin/activate + pip install numpy scipy pytest + PYTHON_PATH=$(which python) + cmake -S . -B "build" -DCMAKE_BUILD_TYPE=Debug -DPYTHON_EXECUTABLE=${PYTHON_PATH} + cd build + make oscode + make -j4 test - Try explicitly including the location of your Eigen library via the - ``CPLUS_INCLUDE_PATH`` environment variable, for example: +Tests for the C++ and python can be run via cmake with .. code:: bash + # This can take a while the first time because of boost + cmake -S . -B "build" -DCMAKE_BUILD_TYPE=Debug -DPYTHON_EXECUTABLE=${PYTHON_PATH} + # If your already inside of the build folder then run + # cmake .. -DCMAKE_BUILD_TYPE=Debug -DPYTHON_EXECUTABLE=${PYTHON_PATH} + cd ./build + make -j4 test - CPLUS_INCLUDE_PATH=/usr/include/eigen3 python setup.py install --user - # or - CPLUS_INCLUDE_PATH=/usr/include/eigen3 pip install pyoscode +Valid `CMAKE_BUILD_TYPE`s are `Release` `Debug`, `LSAN`, `MSAN`, `ASAN`, and `TSAN` + + +Once in the build folder you can run the following to rebuild the cmake if you need to. + + .. code:: bash + cmake .. -DCMAKE_BUILD_TYPE=Debug + make -j4 test - where ``/usr/include/eigen3`` should be replaced with your system-specific - eigen location. Thanks ------ @@ -258,6 +266,8 @@ devs of `exhale `__ for making the beautiful C Changelog --------- +- 1.3.0: + - Add support for cmake build and tests - 1.2.0: - Update the version of Eigen to 3.4.0 - 1.1.2: diff --git a/deps/eigen b/deps/eigen deleted file mode 160000 index e7248b2..0000000 --- a/deps/eigen +++ /dev/null @@ -1 +0,0 @@ -Subproject commit e7248b26a1ed53fa030c5c459f7ea095dfd276ac diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt new file mode 100644 index 0000000..ee71b5e --- /dev/null +++ b/examples/CMakeLists.txt @@ -0,0 +1,25 @@ +cmake_minimum_required(VERSION 3.20.2) +project( + oscode_examples + VERSION 1.3.0 + LANGUAGES C CXX) + + +set(CMAKE_VERBOSE_MAKEFILE YES) +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED YES) +set(CMAKE_CXX_EXTENSIONS NO) + +include(FetchContent) +set(OSCODE_HEADER_ONLY ON) +FetchContent_Declare( + oscode + GIT_REPOSITORY https://github.com/SteveBronder/oscode + GIT_TAG feature/cmake +) + +FetchContent_MakeAvailable(oscode) + +add_executable(burst_ex burst.cpp) + +target_link_libraries(burst_ex PRIVATE oscode) diff --git a/examples/airy.cpp b/examples/airy.cpp index b63b7c3..ff2148e 100644 --- a/examples/airy.cpp +++ b/examples/airy.cpp @@ -153,10 +153,10 @@ int main(){ solution.solve(); /** Extract the solution and the types of steps taken by oscode */ - std::list> xs = solution.sol; - std::list ts = solution.times; - std::list types = solution.wkbs; - std::list> dosol = solution.dosol; + std::vector> xs = solution.sol; + std::vector ts = solution.times; + std::vector types = solution.wkbs; + std::vector> dosol = solution.dosol; int steps = solution.ssteps; /** Write result in file */ diff --git a/examples/burst.cpp b/examples/burst.cpp index e3bd066..c30e33f 100644 --- a/examples/burst.cpp +++ b/examples/burst.cpp @@ -150,9 +150,9 @@ int main(){ solution.solve(); /** Extract the solution and the types of steps taken by oscode */ - std::list> xs = solution.sol; - std::list ts = solution.times; - std::list types = solution.wkbs; + std::vector> xs = solution.sol; + std::vector ts = solution.times; + std::vector types = solution.wkbs; int steps = solution.ssteps; /** Write result in file */ diff --git a/include/oscode/interpolator.hpp b/include/oscode/interpolator.hpp index 11d70d2..359e414 100644 --- a/include/oscode/interpolator.hpp +++ b/include/oscode/interpolator.hpp @@ -2,6 +2,7 @@ #include #include #include +#include template *, typename InputIt_x = double *> diff --git a/include/oscode/rksolver.hpp b/include/oscode/rksolver.hpp index 5d20d2d..a14f1d1 100644 --- a/include/oscode/rksolver.hpp +++ b/include/oscode/rksolver.hpp @@ -3,7 +3,7 @@ #include #include #include -#include +#include /** Class that defines a 4 and 5th order Runge-Kutta method. * @@ -73,9 +73,9 @@ class RKSolver { Eigen::Matrix, 7, 2> k_dense; Eigen::Matrix P_dense; void dense_step(double t0, double h0, std::complex y0, - std::complex dy0, const std::list &dots, - std::list> &doxs, - std::list> &dodxs); + std::complex dy0, const std::vector &dots, + std::vector> &doxs, + std::vector> &dodxs); // Experimental continuous representation of solution Eigen::Matrix, 7, 1> x_vdm; }; @@ -180,9 +180,9 @@ RKSolver::dense_point(std::complex x, std::complex dx, */ void RKSolver::dense_step(double t0, double h0, std::complex y0, std::complex dy0, - const std::list &dots, - std::list> &doxs, - std::list> &dodxs) { + const std::vector &dots, + std::vector> &doxs, + std::vector> &dodxs) { int docount = dots.size(); double h = h0; diff --git a/include/oscode/solver.hpp b/include/oscode/solver.hpp index 4c916da..4a19306 100644 --- a/include/oscode/solver.hpp +++ b/include/oscode/solver.hpp @@ -30,7 +30,7 @@ class Solution { double fend, fnext; /** a boolean encoding the direction of integration: 1/True for forward. */ bool sign; - + bool dense_output_{false}; public: Solution(de_system &de_sys, std::complex x0, std::complex dx0, double t_i, double t_f, int o = 3, double r_tol = 1e-4, @@ -52,24 +52,24 @@ class Solution { int ssteps, totsteps, wkbsteps; /** Lists to contain the solution and its derivative evaluated at internal * points taken by the solver (i.e. not dense output) after a run */ - std::list> sol, dsol; + std::vector> sol, dsol; /** List to contain the timepoints at which the solution and derivative are * internally evaluated by the solver */ - std::list times; + std::vector times; /** List to contain the "type" of each step (RK/WKB) taken internally by the * solver after a run */ - std::list wkbs; + std::vector wkbs; /** Lists to contain the timepoints at which dense output was evaluated. This * list will always be sorted in ascending order (with possible duplicates), * regardless of the order the timepoints were specified upon input. */ - std::list dotimes; + std::vector dotimes; /** Lists to contain the dense output of the solution and its derivative */ - std::list> dosol, dodsol; + std::vector> dosol, dodsol; /** Iterator to iterate over the dense output timepoints, for when these * need to be written out to file */ - std::list::iterator dotit; + //std::vector::iterator dotit; // Experimental: list to contain continuous representation of the solution - std::list, 7, 1>> sol_vdm; + std::vector, 7, 1>> sol_vdm; }; /** Constructor for when dense output was not requested. Sets up solution of the @@ -142,7 +142,7 @@ Solution::Solution(de_system &de_sys, std::complex x0, dotimes.push_back(t_f); dosol.push_back(x0); dodsol.push_back(dx0); - dotit = dotimes.end(); + //dotit = dotimes.end(); switch (order) { case 1: @@ -186,6 +186,7 @@ Solution::Solution(de_system &de_sys, std::complex x0, // Make underlying equation system accessible de_sys_ = &de_sys; + dense_output_ = true; // Set parameters for solver x = x0; dx = dx0; @@ -226,27 +227,22 @@ Solution::Solution(de_system &de_sys, std::complex x0, } // Dense output preprocessing: sort and reverse if necessary - int dosize = do_times.size(); - dotimes.resize(dosize); + std::size_t dosize = do_times.size(); dosol.resize(dosize); dodsol.resize(dosize); // Copy dense output points to list - auto doit = do_times.begin(); - for (auto it = dotimes.begin(); it != dotimes.end(); it++) { - *it = *doit; - doit++; - } + dotimes = do_times; // Sort to ensure ascending order - dotimes.sort(); + std::sort(dotimes.begin(), dotimes.end()); // Reverse if necessary if ((de_sys_->is_interpolated == 1 && de_sys_->Winterp.sign_ == 0) || (de_sys_->is_interpolated == 0 && sign == 0)) { - dotimes.reverse(); + std::reverse(dotimes.begin(), dotimes.end()); } - dotit = dotimes.begin(); + //dotit = dotimes.begin(); switch (order) { case 1: wkbsolver1 = WKBSolver1(*de_sys_, order); @@ -287,7 +283,7 @@ void Solution::solve() { double wkbdelta, rkdelta; std::complex xnext, dxnext; bool wkb = false; - Eigen::Index maxindex_wkb, maxindex_rk; + Eigen::Index maxindex_wkb{0}, maxindex_rk; h = h0; tnext = t + h; // Initialise stats @@ -299,15 +295,14 @@ void Solution::solve() { totsteps = 0; wkbsteps = 0; // Dense output - std::list inner_dotimes; - std::list> inner_dosols, inner_dodsols; - auto it_dosol = dosol.begin(); - auto it_dodsol = dodsol.begin(); + std::vector inner_dotimes; + std::vector> inner_dosols, inner_dodsols; Eigen::Matrix, 1, 2> y_dense_rk; std::complex x_dense_rk, dx_dense_rk; // Experimental continuous solution, vandermonde representation Eigen::Matrix, 7, 1> xvdm; - + std::size_t do_solve_count = 0; + std::size_t do_dodsol_count = 0; while (fend > 0) { // Check if we are reaching the end of integration if (fnext < 0) { @@ -394,13 +389,17 @@ void Solution::solve() { // check if chosen step was successful if (std::abs(hnext) >= std::abs(h)) { // std::cout << "All dense output points: " << std::endl; - if (dotit != dotimes.end()) { + if (dense_output_ && do_solve_count < dotimes.size()) { // std::cout << *dotit << std::endl; - while ((*dotit - t >= 0 && tnext - *dotit >= 0) || - (*dotit - t <= 0 && tnext - *dotit <= 0)) { - inner_dotimes.push_back(*dotit); - dotit++; - } + auto dot_it = dotimes.begin(); + std::advance(dot_it, do_solve_count); + for (; dot_it != dotimes.end() && + ((*dot_it - t >= 0 && tnext - *dot_it >= 0) || + (*dot_it - t <= 0 && tnext - *dot_it <= 0)); + ++dot_it) { + inner_dotimes.push_back(*dot_it); + do_solve_count++; + } if (inner_dotimes.size() > 0) { inner_dosols.resize(inner_dotimes.size()); inner_dodsols.resize(inner_dotimes.size()); @@ -423,16 +422,18 @@ void Solution::solve() { } } } - auto inner_it = inner_dosols.begin(); - auto inner_dit = inner_dodsols.begin(); - while (inner_it != inner_dosols.end() && it_dosol != dosol.end() && - inner_dit != inner_dodsols.end() && it_dodsol != dodsol.end()) { - *it_dosol = *inner_it; - *it_dodsol = *inner_dit; - it_dodsol++; - it_dosol++; - inner_it++; - inner_dit++; + auto it_dosol = dosol.begin(); + std::advance(it_dosol, do_dodsol_count); + auto it_dodsol = dodsol.begin(); + std::advance(it_dodsol, do_dodsol_count); + for (auto inner_it = inner_dosols.begin(), + inner_dit = inner_dodsols.begin(); + inner_it != inner_dosols.end() && it_dosol != dosol.end() && + inner_dit != inner_dodsols.end() && it_dodsol != dodsol.end(); + it_dodsol++, it_dosol++, inner_it++, inner_dit++, + do_dodsol_count++) { + *it_dosol = std::move(*inner_it); + *it_dodsol = std::move(*inner_dit); } inner_dotimes.resize(0); inner_dosols.resize(0); @@ -497,15 +498,15 @@ void Solution::solve() { // reversed at the start) if (de_sys_->is_interpolated == 1) { if (de_sys_->Winterp.sign_ == 0) { - dotimes.reverse(); - dosol.reverse(); - dodsol.reverse(); + std::reverse(dotimes.begin(), dotimes.end()); + std::reverse(dosol.begin(), dosol.end()); + std::reverse(dodsol.begin(), dodsol.end()); } } else { if (sign == 0) { - dotimes.reverse(); - dosol.reverse(); - dodsol.reverse(); + std::reverse(dotimes.begin(), dotimes.end()); + std::reverse(dosol.begin(), dosol.end()); + std::reverse(dodsol.begin(), dodsol.end()); } } diff --git a/include/oscode/system.hpp b/include/oscode/system.hpp index 2d48233..a4ce86c 100644 --- a/include/oscode/system.hpp +++ b/include/oscode/system.hpp @@ -1,6 +1,6 @@ #pragma once #include - +#include /** */ class de_system { private: diff --git a/include/oscode/wkbsolver.hpp b/include/oscode/wkbsolver.hpp index d6b4671..efbac10 100644 --- a/include/oscode/wkbsolver.hpp +++ b/include/oscode/wkbsolver.hpp @@ -3,7 +3,7 @@ #include #include #include -#include +#include /** Class to carry out WKB steps of varying orders. */ class WKBSolver { @@ -103,7 +103,7 @@ class WKBSolver { // error estimate on step std::complex err_fp, err_fm, err_dfp, err_dfm; // dense output - std::list> doxs, dodxs, dows; + std::vector> doxs, dodxs, dows; Eigen::Matrix, 1, 4> dense_s_, dense_ds_, dense_ds_i; std::complex dense_ap_, dense_am_, dense_bp_, dense_bm_; // experimental dense output + quadrature @@ -120,9 +120,9 @@ class WKBSolver { const Eigen::Matrix, 5, 1> &ws5, const Eigen::Matrix, 5, 1> &gs5); // dense output - void dense_step(double t0, const std::list &dots, - std::list> &doxs, - std::list> &dodxs); + void dense_step(double t0, const std::vector &dots, + std::vector> &doxs, + std::vector> &dodxs); Eigen::Matrix dense_weights_6(double t); Eigen::Matrix dense_weights_derivs_6(double t); std::complex @@ -205,8 +205,7 @@ WKBSolver::WKBSolver(de_system &de_sys, int order) { d1w4_5_w << -0.518019493788065, 1.52752523062733, -3.49148624058568, -0.560400043118500e-8, 2.48198050935041; // Set polynomial Gauss--Lobatto coefficients for dense output + quadrature - double a = std::sqrt(147 + 42 * std::sqrt(7.)); - double b = std::sqrt(147 - 42 * std::sqrt(7.)); + interp_vandermonde << 0.06250000000, -0.1946486424, 0.6321486424, 0.6321486424, -0.1946486424, 0.06250000000, -0.06250000000, 0.2544242701, @@ -270,7 +269,6 @@ WKBSolver::step(std::complex x0, std::complex dx0, double t0, Eigen::Matrix, 6, 1> integrand6, s3_interp, s2_interp, s1_interp; Eigen::Matrix, 6, 1> ds1_interp, ds2_interp, ds3_interp; - double t_trans; std::complex s0, s1, s2, s3, dense_fp, dense_fm, dense_x; std::complex ds0, ds1, ds2, ds3, dense_dfpf, dense_dfmf, dense_dx; @@ -427,9 +425,9 @@ WKBSolver::step(std::complex x0, std::complex dx0, double t0, * @param dodxs[in,out] dense output for the derivative of the solution * \f$\dot{x}\f$ */ -void WKBSolver::dense_step(double t0, const std::list &dots, - std::list> &doxs, - std::list> &dodxs) { +void WKBSolver::dense_step(double t0, const std::vector &dots, + std::vector> &doxs, + std::vector> &dodxs) { // We have: ws_, gs_, ws5_, gs5_, ws7_, x, dx, ddx, h, dws_, dws5_, d2wx, // d3wx, etc., diff --git a/pyoscode/_pyoscode.cpp b/pyoscode/_pyoscode.cpp index c13b43d..a1a5216 100644 --- a/pyoscode/_pyoscode.cpp +++ b/pyoscode/_pyoscode.cpp @@ -104,7 +104,7 @@ static PyObject *_pyoscode_solve_fn(PyObject *self, PyObject *args, PyObject *kw t_eval = (double*)PyArray_DATA(t_evalarray_arr); t_evalsize = (int)PyArray_SIZE(t_evalarray_arr); } - std::list t_evallist; + std::vector t_evallist; t_evallist.resize(t_evalsize); int i=0; for(auto it=t_evallist.begin(); it!=t_evallist.end(); it++){ @@ -134,11 +134,11 @@ static PyObject *_pyoscode_solve_fn(PyObject *self, PyObject *args, PyObject *kw return (PyObject *) NULL; } de_system sys = de_system(&wfun, &gfun); - std::list> sol,dsol; - std::list times; - std::list wkbs; - std::list> x_eval, dx_eval; - std::list,7,1>> cts_rep; + std::vector> sol,dsol; + std::vector times; + std::vector wkbs; + std::vector> x_eval, dx_eval; + std::vector,7,1>> cts_rep; if(t_evalobj!=NULL){ Solution solution(sys,x0,dx0,ti,tf,t_evallist,order,rtol,atol,h0,full_output); @@ -269,7 +269,7 @@ static PyObject *_pyoscode_solve(PyObject *self, PyObject *args, PyObject *kwarg t_eval = (double*)PyArray_DATA(t_evalarray_arr); t_evalsize = (int)PyArray_SIZE(t_evalarray_arr); } - std::list t_evallist; + std::vector t_evallist; t_evallist.resize(t_evalsize); int i=0; for(auto it=t_evallist.begin(); it!=t_evallist.end(); it++){ @@ -346,11 +346,11 @@ regions of the solution efficiently and numerical inaccuracies. Please \ consider refining the sampling of the array(s) or switching to a more \ suitable independent variable.",1); } - std::list> sol,dsol; - std::list times; - std::list wkbs; - std::list> x_eval, dx_eval; - std::list,7,1>> cts_rep; + std::vector> sol,dsol; + std::vector times; + std::vector wkbs; + std::vector> x_eval, dx_eval; + std::vector,7,1>> cts_rep; if(t_evalobj!=NULL){ Solution solution(sys,x0,dx0,ti,tf,t_evallist,order,rtol,atol,h0,full_output); diff --git a/setup.py b/setup.py index 0e3afd1..1da7816 100644 --- a/setup.py +++ b/setup.py @@ -1,8 +1,11 @@ from __future__ import absolute_import, with_statement, print_function, division from setuptools import setup, Extension, find_packages +import sys import os +import platform import numpy as np +source_dir = os.getenv('OSCODE_EIGEN_INCLUDE_DIR') def readme(short=False): with open("README.rst") as f: if short: @@ -10,17 +13,28 @@ def readme(short=False): else: return f.read() +extra_compile_args = [] +if platform.system() == 'Windows': # For Windows + if "MSC" in platform.python_compiler(): + # Visual Studio (MSVC) + extra_compile_args = ['/std:c++17'] + else: + # assume MinGW or similar + extra_compile_args = ['-std=c++17', '-Wall'] +else: # For Unix/Linux/MacOS + extra_compile_args = ['-std=c++17', '-Wall'] + pyoscode_module = Extension( name="_pyoscode", sources=["pyoscode/_pyoscode.cpp"], - include_dirs=['include','pyoscode',np.get_include(),'deps/eigen'], + include_dirs=['include','pyoscode',np.get_include(), source_dir], depends=["pyoscode/_python.hpp", "pyoscode/_pyoscode.hpp"], - extra_compile_args=['-std=c++17','-Wall'] + extra_compile_args=extra_compile_args ) setup( name="pyoscode", - version="1.2.0", + version="1.3.0", description=readme(short=True), long_description=readme(), url="https://github.com/fruzsinaagocs/oscode", diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt new file mode 100644 index 0000000..e7503a4 --- /dev/null +++ b/tests/CMakeLists.txt @@ -0,0 +1,22 @@ +add_executable(burst_test burst_test.cpp) +target_include_directories(burst_test PRIVATE ${oscode_INCLUDE_DIR}) +target_link_libraries(burst_test PRIVATE Eigen3::Eigen gtest gtest_main) + +add_executable(aries_test aries_test.cpp) +target_include_directories(aries_test PRIVATE ${oscode_INCLUDE_DIR} ${boost_INCLUDE_DIR}) +target_link_libraries(aries_test PRIVATE Eigen3::Eigen gtest gtest_main) + +enable_testing() +if (WIN32 AND NOT CMAKE_CONFIGURATION_TYPES) + set(config_subdir "/$") +else() + set(config_subdir "") +endif() + +add_test(NAME burst COMMAND ${CMAKE_CURRENT_BINARY_DIR}${config_subdir}/burst_test) + +find_package(Python COMPONENTS Interpreter REQUIRED) +add_test(NAME python_tests COMMAND ${PYTHON_EXECUTABLE} -m pytest ${CMAKE_CURRENT_SOURCE_DIR}) + +add_custom_target(test COMMAND ${CMAKE_CTEST_COMMAND} -C ${CMAKE_BUILD_TYPE} --extra-verbose DEPENDS burst_test pyoscode) + diff --git a/tests/aries_test.cpp b/tests/aries_test.cpp new file mode 100644 index 0000000..5a39dc4 --- /dev/null +++ b/tests/aries_test.cpp @@ -0,0 +1,176 @@ +#include +#include +#include +#include +#include +#include +#include + +/** Defines the friction term in the ODE to solve + * @param[in] t Value of the independent variable in the ODE + * @returns The value of the friction term at \a t + */ +std::complex g(double t) { return 0.0; } + +/** Defines the frequency term in the ODE to solve + * @param[in] t Value of the independent variable in the ODE + * @returns The value of the frequency term at \a t + */ +std::complex w(double t) { return std::sqrt(t); } + +/** Defines the initial value of the solution of the ODE + * @param[in] t Value of the independent variable in the ODE + * @returns The value of the solution \a x at \a t + */ +std::complex xairy(double t) { + return std::complex(boost::math::airy_ai(-t), + boost::math::airy_bi(-t)); +} + +/** Defines the initial value of the derivative of the solution of the ODE + * @param[in] t Value of the independent variable in the ODE + * @returns The value of the derivative of solution \a dx/dt at \a t + */ +std::complex dxairy(double t) { + return std::complex(-boost::math::airy_ai_prime(-t), + -boost::math::airy_bi_prime(-t)); +} + +/** \brief Routine to call oscode to solve an example ODE, to illustrate basic + * functionality. + * + * Routine to call oscode to solve the example ODE (the "burst" equation): + * \f[ + * \ddot{x} + \fract{n^2-1}{(1+t^2)^2}x = 0 + * \f] + * where \f$ n \f$ is a parameter that controls how many oscillations the + * solution \f$ x(t) \f$ will go through. + * + * The routine defines the differential equation in four different ways: + * -# Defines the frequency and friction terms (\f$ \omega(t) \f$ and \f$ + * \gamma(t) \f$) as functions, + * -# Defines the frequency and friction terms as sequences: + * -# as Eigen vectors, + * -# as std::arrays, + * -# as std::vectors, + * and then feeds a pointer to the underlying data array to the class \a + * de_system representing the ODE. + * + * All four methods are included in this routine, with the last three commented + * out. Test any of the above methods by uncommenting and commenting out all + * others. + * + * After solving the ODE, this routine prints the solution and its derivative to + * a file called \a output.txt, the contents of which you can plot with the + * Python script \a plot_burst.py, or any other script of your choice. + */ +TEST(SolverTest, SolveAriesUnEvenFwd) { + /** Define integration range */ + double ti = 1.0; + double tf = 1e6; + + /** Define initial conditions */ + std::complex x0 = xairy(ti); + std::complex dx0 = dxairy(ti); + + /** We'll ask for dense output at the following points */ + + double a = 0.1; + double t_dense_i = 5e1; + std::vector t_dense(200); + std::generate( + t_dense.begin(), t_dense.end(), + [n = 0, &a, &t_dense_i]() mutable { return n++ * a + t_dense_i; }); + + /** Create differential equation "system" */ + /** Method 1: Give frequency and damping term as functions */ + de_system sys(&w, &g); + /** Solve the ODE */ + Solution solution(sys, x0, dx0, ti, tf, t_dense, 3, 1e-8); + solution.solve(); + std::vector> sol_vec; + sol_vec.reserve(t_dense.size()); + for (auto &sol : solution.dosol) { + sol_vec.push_back(sol); + } + std::vector> true_sol_vec; + true_sol_vec.reserve(t_dense.size()); + for (auto &t : t_dense) { + true_sol_vec.push_back(xairy(t)); + } + std::vector> dsol_vec; + for (auto &dsol : solution.dodsol) { + dsol_vec.push_back(dsol); + } + std::vector> true_dsol_vec; + true_dsol_vec.reserve(t_dense.size()); + for (auto &t : t_dense) { + true_dsol_vec.push_back(dxairy(t)); + } + for (std::size_t i = 0; i < t_dense.size(); ++i) { + EXPECT_NEAR((std::real(sol_vec[i]) - std::real(true_sol_vec[i])), 0.0f, + 1e-4f); + EXPECT_NEAR((std::imag(sol_vec[i]) - std::imag(true_sol_vec[i])), 0.0f, + 1e-4f); + EXPECT_NEAR((std::real(dsol_vec[i]) - std::real(true_dsol_vec[i])), 0.0f, + 1e-4f); + EXPECT_NEAR((std::imag(dsol_vec[i]) - std::imag(true_dsol_vec[i])), 0.0f, + 1e-4f); + + } +} +/* +TEST(SolverTest, SolveAriesUnEvenBwd) { + double ti = 1.0; + double tf = 1e6; + std::complex x0 = xairy(ti); + std::complex dx0 = dxairy(ti); + + double a = 0.1; + double t_dense_i = 5e1; + std::vector t_dense(200); + std::generate( + t_dense.begin(), t_dense.end(), + [n = 0, &a, &t_dense_i]() mutable { return n++ * a + t_dense_i; }); + std::reverse(t_dense.begin(), t_dense.end()); + + de_system sys(&w, &g); + Solution solution(sys, x0, dx0, tf, ti, t_dense, 3, 1e-8, 0, -0.1); + solution.solve(); + std::vector> sol_vec; + for (auto &sol : solution.dosol) { + sol_vec.push_back(sol); + } + std::vector> true_sol_vec; + for (auto &t : t_dense) { + true_sol_vec.push_back(xairy(t)); + } + std::vector> dsol_vec; + for (auto &dsol : solution.dodsol) { + dsol_vec.push_back(dsol); + } + + std::vector> true_dsol_vec; + for (auto &t : t_dense) { + true_dsol_vec.push_back(dxairy(t)); + } + for (std::size_t i = 0; i < t_dense.size(); ++i) { + EXPECT_NEAR((std::real(sol_vec[i]) - std::real(true_sol_vec[i])), 0.0f, + 1e-3f) + << "i = " << i << " sol = " << sol_vec[i] + << " true_sol = " << true_sol_vec[i]; + EXPECT_NEAR((std::imag(sol_vec[i]) + std::imag(true_sol_vec[i])), 0.0f, + 1e-4f) + << "i = " << i << " sol = " << sol_vec[i] + << " true_sol = " << true_sol_vec[i]; + EXPECT_NEAR((std::real(dsol_vec[i]) + std::real(true_dsol_vec[i])), 0.0f, + 1e-4f) + << "i = " << i << " dsol = " << dsol_vec[i] + << " true_dsol = " << true_dsol_vec[i]; + EXPECT_NEAR((std::imag(dsol_vec[i]) - std::imag(true_dsol_vec[i])), 0.0f, + 1e-4f) + << "i = " << i << " dsol = " << dsol_vec[i] + << " true_dsol = " << true_dsol_vec[i]; + } +} +*/ \ No newline at end of file diff --git a/tests/burst_test.cpp b/tests/burst_test.cpp new file mode 100644 index 0000000..ccf9b89 --- /dev/null +++ b/tests/burst_test.cpp @@ -0,0 +1,191 @@ +#include +#include +#include +#include +#include +#include + +/** A variable to control the number of oscillations in the solution of the + * burst equation */ +double n = 100.0; + +/** Defines the friction term in the ODE to solve + * @param[in] t Value of the independent variable in the ODE + * @returns The value of the friction term at \a t + */ +std::complex g(double t) { return 0.0; }; + +/** Defines the frequency term in the ODE to solve + * @param[in] t Value of the independent variable in the ODE + * @returns The value of the frequency term at \a t + */ +std::complex w(double t) { + return std::pow(n * n - 1.0, 0.5) / (1.0 + t * t); +}; + +/** Defines the initial value of the solution of the ODE + * @param[in] t Value of the independent variable in the ODE + * @returns The value of the solution \a x at \a t + */ +std::complex xburst(double t) { + return 100 * std::pow(1.0 + t * t, 0.5) / n * + std::complex(std::cos(n * std::atan(t)), + std::sin(n * std::atan(t))); +}; + +/** Defines the initial value of the derivative of the solution of the ODE + * @param[in] t Value of the independent variable in the ODE + * @returns The value of the derivative of solution \a dx/dt at \a t + */ +std::complex dxburst(double t) { + return 100 / std::pow(1.0 + t * t, 0.5) / n * + (std::complex(t, n) * std::cos(n * std::atan(t)) + + std::complex(-n, t) * std::sin(n * std::atan(t))); +}; + +template +std::vector linspace(T start, T end, std::size_t points) { + std::vector res(points); + float step = (end - start) / (points - 1); + size_t i = 0; + for (auto &e : res) { + e = start + step * i++; + } + return res; +} + +/** \brief Routine to call oscode to solve an example ODE, to illustrate basic + * functionality. + * + * Routine to call oscode to solve the example ODE (the "burst" equation): + * \f[ + * \ddot{x} + \fract{n^2-1}{(1+t^2)^2}x = 0 + * \f] + * where \f$ n \f$ is a parameter that controls how many oscillations the + * solution \f$ x(t) \f$ will go through. + * + * The routine defines the differential equation in four different ways: + * -# Defines the frequency and friction terms (\f$ \omega(t) \f$ and \f$ + * \gamma(t) \f$) as functions, + * -# Defines the frequency and friction terms as sequences: + * -# as Eigen vectors, + * -# as std::arrays, + * -# as std::vectors, + * and then feeds a pointer to the underlying data array to the class \a + * de_system representing the ODE. + * + * All four methods are included in this routine, with the last three commented + * out. Test any of the above methods by uncommenting and commenting out all + * others. + * + * After solving the ODE, this routine prints the solution and its derivative to + * a file called \a output.txt, the contents of which you can plot with the + * Python script \a plot_burst.py, or any other script of your choice. + */ + +TEST(SolverTest, SolveBurstEvenFwd) { + /** Define integration range */ + double ti = -2 * n; + double tf = 2 * n; + + /** Define initial conditions */ + std::complex x0 = xburst(ti); + std::complex dx0 = dxburst(ti); + + /** Create differential equation "system" */ + /** Method 1: Give frequency and damping term as functions */ + de_system sys(&w, &g); + /** Solve the ODE */ + std::vector times = linspace(ti, tf, 1000); + Solution solution(sys, x0, dx0, ti, tf, times, 3, 1e-8); + solution.solve(); + std::vector time_comp_vec = linspace(ti, tf, 1000); + std::vector> sol_vec; + for (auto &sol : solution.dosol) { + sol_vec.push_back(sol); + } + std::vector> true_sol_vec; + for (auto &t : time_comp_vec) { + true_sol_vec.push_back(xburst(t)); + } + std::vector> dsol_vec; + for (auto &dsol : solution.dodsol) { + dsol_vec.push_back(dsol); + } + + std::vector> true_dsol_vec; + for (auto &t : time_comp_vec) { + true_dsol_vec.push_back(dxburst(t)); + } + for (int i = 0; i < time_comp_vec.size(); ++i) { + EXPECT_NEAR((std::real(sol_vec[i]) - std::real(true_sol_vec[i])), 0.0f, + 1e-2f); + EXPECT_NEAR((std::imag(sol_vec[i]) - std::imag(true_sol_vec[i])), 0.0f, + 1e-2f); + EXPECT_NEAR((std::real(dsol_vec[i]) - std::real(true_dsol_vec[i])), 0.0f, + 1e-2); + EXPECT_NEAR((std::imag(dsol_vec[i]) - std::imag(true_dsol_vec[i])), 0.0f, + 1e-2); + // std::cout << true_sol_vec[i] << ", " << sol_vec[i] << ": " << + // sol_vec[i] - true_sol_vec[i] << std::endl; std::cout << + // true_dsol_vec[i] << ", " << dsol_vec[i] << ": " << dsol_vec[i] - + // true_dsol_vec[i] << std::endl; + } +} + +TEST(SolverTest, SolveBurstEvenRev) { + /** Define integration range */ + double ti = -2 * n; + double tf = 2 * n; + + /** Define initial conditions */ + std::complex x0 = xburst(tf); + std::complex dx0 = dxburst(tf); + + /** Create differential equation "system" */ + /** Method 1: Give frequency and damping term as functions */ + de_system sys(&w, &g); + /** Solve the ODE */ + std::vector times = linspace(ti, tf, 1000); + std::reverse(times.begin(), times.end()); + + Solution solution(sys, x0, dx0, tf, ti, times, 3, 1e-12, 0, -1); + solution.solve(); + std::vector> sol_vec; + for (auto &sol : solution.dosol) { + sol_vec.push_back(sol); + } + std::vector time_comp_vec = linspace(ti, tf, 1000); + std::reverse(time_comp_vec.begin(), time_comp_vec.end()); + std::vector> true_sol_vec; + for (auto &t : time_comp_vec) { + true_sol_vec.push_back(xburst(t)); + } + std::vector> dsol_vec; + for (auto &dsol : solution.dodsol) { + dsol_vec.push_back(dsol); + } + + std::vector> true_dsol_vec; + for (auto &t : time_comp_vec) { + true_dsol_vec.push_back(dxburst(t)); + } + for (int i = 0; i < time_comp_vec.size(); ++i) { + EXPECT_NEAR((std::real(sol_vec[i]) - std::real(true_sol_vec[i])), 0.0f, + 1e-3f) + << "i = " << i << " sol = " << sol_vec[i] + << " true_sol = " << true_sol_vec[i]; + EXPECT_NEAR((std::imag(sol_vec[i]) + std::imag(true_sol_vec[i])), 0.0f, + 1e-4f) + << "i = " << i << " sol = " << sol_vec[i] + << " true_sol = " << true_sol_vec[i]; + EXPECT_NEAR((std::real(dsol_vec[i]) + std::real(true_dsol_vec[i])), 0.0f, + 1e-4f) + << "i = " << i << " dsol = " << dsol_vec[i] + << " true_dsol = " << true_dsol_vec[i]; + EXPECT_NEAR((std::imag(dsol_vec[i]) - std::imag(true_dsol_vec[i])), 0.0f, + 1e-4f) + << "i = " << i << " dsol = " << dsol_vec[i] + << " true_dsol = " << true_dsol_vec[i]; + } +} diff --git a/tests/test_airy.py b/tests/test_airy.py index 141f517..44240ee 100644 --- a/tests/test_airy.py +++ b/tests/test_airy.py @@ -21,20 +21,20 @@ def test_no_integration(): def test_airy_forward_even(): # Integration forward on even grid - ts = np.linspace(1,100,5000) - ws = np.sqrt(ts) - gs = np.zeros_like(ws) - ti = 1.0 - tf = 100.0 - x0 = airy(-ti)[0] + 1j*airy(-ti)[2] - dx0 = -airy(-ti)[1] - 1j*airy(-ti)[3] - t_eval = np.linspace(ti,tf,5000) - sol = pyoscode.solve(ts, ws, gs, ti, tf, x0, dx0, t_eval=t_eval, even_grid=True) + time_space = np.linspace(1,100,5000) + weights = np.sqrt(time_space) + gradients = np.zeros_like(weights) + initial_timepoint = 1.0 + final_timepoint = 100.0 + initial_value = airy(-initial_timepoint)[0] + 1j*airy(-initial_timepoint)[2] + initial_deriv_value = -airy(-initial_timepoint)[1] - 1j*airy(-initial_timepoint)[3] + t_eval = np.linspace(initial_timepoint,final_timepoint,5000) + sol = pyoscode.solve(time_space, weights, gradients, initial_timepoint, final_timepoint, initial_value, initial_deriv_value, t_eval=t_eval, even_grid=True) t = np.asarray(sol['t']) x = np.asarray(sol['sol']) dense = np.asarray(sol['x_eval']) dense_d = np.asarray(sol['dx_eval']) - ana_t = np.linspace(ti,tf,5000) + ana_t = np.linspace(initial_timepoint,final_timepoint,5000) ana_x = np.asarray([airy(-T)[0]+1j*airy(-T)[2] for T in t]) dense_ana_x = np.asarray([airy(-T)[0]+1j*airy(-T)[2] for T in ana_t]) dense_ana_dx = np.asarray([-airy(-T)[1]-1j*airy(-T)[3] for T in ana_t]) diff --git a/tests/test_arrays.py b/tests/test_arrays.py index bdc46d5..0142713 100644 --- a/tests/test_arrays.py +++ b/tests/test_arrays.py @@ -97,6 +97,4 @@ def test_not_monotonous(): dx0 = 0.0 with pytest.raises(TypeError): pyoscode.solve(ts,ws,gs,ti,tf,x0,dx0) - return None; - - + return None; \ No newline at end of file