From cc625946786b8e2fa39d0585738b82dc884207b5 Mon Sep 17 00:00:00 2001 From: mattoverby Date: Wed, 8 Oct 2025 14:08:24 -0500 Subject: [PATCH 1/2] update lbfgs and cmake --- CMakeLists.txt | 58 +-- LICENSE.txt | 2 +- cmake/libigl.cmake | 11 + deps/mcloptlib/CMakeLists.txt | 42 -- deps/mcloptlib/LICENSE | 21 - deps/mcloptlib/README.md | 29 -- deps/mcloptlib/cmake/FindEigen3.cmake | 81 ---- deps/mcloptlib/include/MCL/Backtracking.hpp | 150 ------ deps/mcloptlib/include/MCL/LBFGS.hpp | 158 ------- deps/mcloptlib/include/MCL/Minimizer.hpp | 113 ----- deps/mcloptlib/include/MCL/MoreThuente.hpp | 330 ------------- deps/mcloptlib/include/MCL/Newton.hpp | 79 ---- deps/mcloptlib/include/MCL/NonLinearCG.hpp | 86 ---- deps/mcloptlib/include/MCL/Problem.hpp | 135 ------ deps/mcloptlib/include/MCL/TrustRegion.hpp | 208 --------- deps/mcloptlib/include/MCL/WolfeBisection.hpp | 110 ----- deps/mcloptlib/test/TestProblem.hpp | 130 ------ deps/mcloptlib/test/testSolvers.cpp | 210 --------- src/Interp.hpp | 2 +- src/LBFGS.hpp | 440 ++++++++++++++++++ src/MPM.cpp | 1 + src/Particle.hpp | 1 + src/Solver.cpp | 47 +- src/Solver.hpp | 9 +- test/solver.cpp | 9 + 25 files changed, 507 insertions(+), 1955 deletions(-) create mode 100644 cmake/libigl.cmake delete mode 100644 deps/mcloptlib/CMakeLists.txt delete mode 100644 deps/mcloptlib/LICENSE delete mode 100644 deps/mcloptlib/README.md delete mode 100644 deps/mcloptlib/cmake/FindEigen3.cmake delete mode 100644 deps/mcloptlib/include/MCL/Backtracking.hpp delete mode 100644 deps/mcloptlib/include/MCL/LBFGS.hpp delete mode 100644 deps/mcloptlib/include/MCL/Minimizer.hpp delete mode 100644 deps/mcloptlib/include/MCL/MoreThuente.hpp delete mode 100644 deps/mcloptlib/include/MCL/Newton.hpp delete mode 100644 deps/mcloptlib/include/MCL/NonLinearCG.hpp delete mode 100644 deps/mcloptlib/include/MCL/Problem.hpp delete mode 100644 deps/mcloptlib/include/MCL/TrustRegion.hpp delete mode 100644 deps/mcloptlib/include/MCL/WolfeBisection.hpp delete mode 100644 deps/mcloptlib/test/TestProblem.hpp delete mode 100644 deps/mcloptlib/test/testSolvers.cpp create mode 100644 src/LBFGS.hpp create mode 100644 test/solver.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 7f925de..8fc6554 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -25,6 +25,7 @@ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") set(CMAKE_BUILD_TYPE Release) add_definitions(-DMPM_SRC_DIR="${CMAKE_CURRENT_SOURCE_DIR}") +option(MPM_TESTING_ONLY "Only compile lib and tests" OFF) ############################################################ # @@ -32,25 +33,8 @@ add_definitions(-DMPM_SRC_DIR="${CMAKE_CURRENT_SOURCE_DIR}") # ############################################################ -# libigl for rendering -include(DownloadProject) -option(LIBIGL_USE_STATIC_LIBRARY "Use libigl as static library" OFF) # off = header only -option(LIBIGL_WITH_COMISO "Use CoMiso" OFF) -option(LIBIGL_WITH_EMBREE "Use Embree" OFF) -option(LIBIGL_WITH_OPENGL "Use OpenGL" ON) -option(LIBIGL_WITH_OPENGL_GLFW "Use GLFW" ON) -option(LIBIGL_WITH_OPENGL_GLFW_IMGUI "Use ImGui" OFF) -option(LIBIGL_WITH_PNG "Use PNG" OFF) -option(LIBIGL_WITH_TETGEN "Use Tetgen" OFF) -option(LIBIGL_WITH_TRIANGLE "Use Triangle" OFF) -option(LIBIGL_WITH_PREDICATES "Use exact predicates" OFF) -option(LIBIGL_WITH_XML "Use XML" OFF) -download_project(PROJ libigl - GIT_REPOSITORY https://github.com/libigl/libigl.git - GIT_TAG main - UPDATE_DISCONNECTED 1 - QUIET) -add_subdirectory(${libigl_SOURCE_DIR} ${libigl_BINARY_DIR}) +# Libigl also includes Eigen, etc +include(libigl) # OpenMP find_package(OpenMP) @@ -60,20 +44,18 @@ if (OPENMP_FOUND) add_definitions(-DOMP_NESTED) endif() -# Include headers -include_directories(${CMAKE_CURRENT_SOURCE_DIR}/src) -include_directories(SYSTEM ${CMAKE_CURRENT_SOURCE_DIR}/deps/mcloptlib/include) -include_directories(SYSTEM ${glad_SOURCE_DIR}/include) -include_directories(SYSTEM ${libigl_SOURCE_DIR}/include) -include_directories(SYSTEM ${libigl_SOURCE_DIR}/external/eigen) +set(MPM_SRC + ${CMAKE_CURRENT_SOURCE_DIR}/src/LBFGS.hpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/Solver.hpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/Solver.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/MPM.hpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/MPM.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/Interp.hpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/Particle.hpp) -set(MPM_SRCS - src/Solver.hpp - src/Solver.cpp - src/MPM.hpp - src/MPM.cpp - src/Interp.hpp - src/Particle.hpp) + +add_library(mpm_optimization ${MPM_SRC}) +target_link_libraries(mpm_optimization PUBLIC igl::core) ############################################################ # @@ -81,5 +63,13 @@ set(MPM_SRCS # ############################################################ -add_executable(sphere src/sphere.cpp ${MPM_SRCS}) -target_link_libraries(sphere glad glfw) +if (NOT MPM_TESTING_ONLY) + igl_include(glfw) + add_executable(sphere ${CMAKE_CURRENT_SOURCE_DIR}/src/sphere.cpp) + target_link_libraries(sphere PUBLIC mpm_optimization igl::glfw) +endif() + +enable_testing() +add_executable(test_solver ${CMAKE_CURRENT_SOURCE_DIR}/test/solver.cpp) +target_link_libraries(test_solver PUBLIC mpm_optimization) +add_test(NAME TestSolver COMMAND test_solver) \ No newline at end of file diff --git a/LICENSE.txt b/LICENSE.txt index 56ef4ca..13de209 100644 --- a/LICENSE.txt +++ b/LICENSE.txt @@ -10,7 +10,7 @@ of conditions and the following disclaimer in the documentation and/or other mat provided with the distribution. THIS SOFTWARE IS PROVIDED "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF MINNESOTA, DULUTH OR CONTRIBUTORS BE +ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF MINNESOTA OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER diff --git a/cmake/libigl.cmake b/cmake/libigl.cmake new file mode 100644 index 0000000..1aae978 --- /dev/null +++ b/cmake/libigl.cmake @@ -0,0 +1,11 @@ +if(TARGET igl::core) + return() +endif() + +include(FetchContent) +FetchContent_Declare( + libigl + GIT_REPOSITORY https://github.com/libigl/libigl.git + GIT_TAG v2.5.0 +) +FetchContent_MakeAvailable(libigl) \ No newline at end of file diff --git a/deps/mcloptlib/CMakeLists.txt b/deps/mcloptlib/CMakeLists.txt deleted file mode 100644 index d33d831..0000000 --- a/deps/mcloptlib/CMakeLists.txt +++ /dev/null @@ -1,42 +0,0 @@ -# The MIT License (MIT) -# Copyright (c) 2017 Matt Overby -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. - -cmake_minimum_required(VERSION 3.1) -project(mcloptlib C CXX) -set(CMAKE_CXX_STANDARD 11) -set(CMAKE_CXX_STANDARD_REQUIRED ON) -set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH}) -set(CMAKE_BUILD_TYPE Debug) -add_definitions( -DMCL_DEBUG=1 ) - -if(CMAKE_COMPILER_IS_GNUCC OR CMAKE_COMPILER_IS_GNUCXX) - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Wno-long-long") -endif() - -find_package(Eigen3 REQUIRED) -include_directories(SYSTEM ${EIGEN3_INCLUDE_DIR}) -include_directories(${CMAKE_CURRENT_SOURCE_DIR}/include) - -enable_testing() -add_executable(testSolvers test/testSolvers.cpp) -add_test(testLBFGS testSolvers lbfgs) -add_test(testCG testSolvers cg) -add_test(testNewton testSolvers newton) diff --git a/deps/mcloptlib/LICENSE b/deps/mcloptlib/LICENSE deleted file mode 100644 index c95244f..0000000 --- a/deps/mcloptlib/LICENSE +++ /dev/null @@ -1,21 +0,0 @@ -The MIT License (MIT) - -Copyright (c) 2017 Matt Overby - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. diff --git a/deps/mcloptlib/README.md b/deps/mcloptlib/README.md deleted file mode 100644 index bc89cbb..0000000 --- a/deps/mcloptlib/README.md +++ /dev/null @@ -1,29 +0,0 @@ -# mcloptlib - -By Matt Overby -[http://www.mattoverby.net](http://www.mattoverby.net) - -mcloptlib is a header-only optimization library for C++ using Eigen and is geared towards lower-dimension graphics problems. -Originally a fork of [Patrick Wieschollek's CppOptimizationLibrary](https://github.com/PatWie/CppNumericalSolvers), but has diverged considerably. - -## Contents: - -Optimization algorithms: -- Newton's -- Non-linear conjugate gradient -- L-BFGS -- Trust Region with - - Cauchy Point - - Dog Leg - -Linesearch methods: -- Backtracking (Armijo) -- Backtracking with cubic interpolation -- Bisection -- MoreThuente - -## To-do: - -- Option of std::function for value/gradient instead of Problem class -- Sparse Hessians -- Auto-diff diff --git a/deps/mcloptlib/cmake/FindEigen3.cmake b/deps/mcloptlib/cmake/FindEigen3.cmake deleted file mode 100644 index 9c546a0..0000000 --- a/deps/mcloptlib/cmake/FindEigen3.cmake +++ /dev/null @@ -1,81 +0,0 @@ -# - Try to find Eigen3 lib -# -# This module supports requiring a minimum version, e.g. you can do -# find_package(Eigen3 3.1.2) -# to require version 3.1.2 or newer of Eigen3. -# -# Once done this will define -# -# EIGEN3_FOUND - system has eigen lib with correct version -# EIGEN3_INCLUDE_DIR - the eigen include directory -# EIGEN3_VERSION - eigen version - -# Copyright (c) 2006, 2007 Montel Laurent, -# Copyright (c) 2008, 2009 Gael Guennebaud, -# Copyright (c) 2009 Benoit Jacob -# Redistribution and use is allowed according to the terms of the 2-clause BSD license. - -if(NOT Eigen3_FIND_VERSION) - if(NOT Eigen3_FIND_VERSION_MAJOR) - set(Eigen3_FIND_VERSION_MAJOR 2) - endif(NOT Eigen3_FIND_VERSION_MAJOR) - if(NOT Eigen3_FIND_VERSION_MINOR) - set(Eigen3_FIND_VERSION_MINOR 91) - endif(NOT Eigen3_FIND_VERSION_MINOR) - if(NOT Eigen3_FIND_VERSION_PATCH) - set(Eigen3_FIND_VERSION_PATCH 0) - endif(NOT Eigen3_FIND_VERSION_PATCH) - - set(Eigen3_FIND_VERSION "${Eigen3_FIND_VERSION_MAJOR}.${Eigen3_FIND_VERSION_MINOR}.${Eigen3_FIND_VERSION_PATCH}") -endif(NOT Eigen3_FIND_VERSION) - -macro(_eigen3_check_version) - file(READ "${EIGEN3_INCLUDE_DIR}/Eigen/src/Core/util/Macros.h" _eigen3_version_header) - - string(REGEX MATCH "define[ \t]+EIGEN_WORLD_VERSION[ \t]+([0-9]+)" _eigen3_world_version_match "${_eigen3_version_header}") - set(EIGEN3_WORLD_VERSION "${CMAKE_MATCH_1}") - string(REGEX MATCH "define[ \t]+EIGEN_MAJOR_VERSION[ \t]+([0-9]+)" _eigen3_major_version_match "${_eigen3_version_header}") - set(EIGEN3_MAJOR_VERSION "${CMAKE_MATCH_1}") - string(REGEX MATCH "define[ \t]+EIGEN_MINOR_VERSION[ \t]+([0-9]+)" _eigen3_minor_version_match "${_eigen3_version_header}") - set(EIGEN3_MINOR_VERSION "${CMAKE_MATCH_1}") - - set(EIGEN3_VERSION ${EIGEN3_WORLD_VERSION}.${EIGEN3_MAJOR_VERSION}.${EIGEN3_MINOR_VERSION}) - if(${EIGEN3_VERSION} VERSION_LESS ${Eigen3_FIND_VERSION}) - set(EIGEN3_VERSION_OK FALSE) - else(${EIGEN3_VERSION} VERSION_LESS ${Eigen3_FIND_VERSION}) - set(EIGEN3_VERSION_OK TRUE) - endif(${EIGEN3_VERSION} VERSION_LESS ${Eigen3_FIND_VERSION}) - - if(NOT EIGEN3_VERSION_OK) - - message(STATUS "Eigen3 version ${EIGEN3_VERSION} found in ${EIGEN3_INCLUDE_DIR}, " - "but at least version ${Eigen3_FIND_VERSION} is required") - endif(NOT EIGEN3_VERSION_OK) -endmacro(_eigen3_check_version) - -if (EIGEN3_INCLUDE_DIR) - - # in cache already - _eigen3_check_version() - set(EIGEN3_FOUND ${EIGEN3_VERSION_OK}) - -else (EIGEN3_INCLUDE_DIR) - - find_path(EIGEN3_INCLUDE_DIR NAMES signature_of_eigen3_matrix_library - PATHS - ${CMAKE_INSTALL_PREFIX}/include - ${KDE4_INCLUDE_DIR} - PATH_SUFFIXES eigen3 eigen - ) - - if(EIGEN3_INCLUDE_DIR) - _eigen3_check_version() - endif(EIGEN3_INCLUDE_DIR) - - include(FindPackageHandleStandardArgs) - find_package_handle_standard_args(Eigen3 DEFAULT_MSG EIGEN3_INCLUDE_DIR EIGEN3_VERSION_OK) - - mark_as_advanced(EIGEN3_INCLUDE_DIR) - -endif(EIGEN3_INCLUDE_DIR) - diff --git a/deps/mcloptlib/include/MCL/Backtracking.hpp b/deps/mcloptlib/include/MCL/Backtracking.hpp deleted file mode 100644 index 90d0764..0000000 --- a/deps/mcloptlib/include/MCL/Backtracking.hpp +++ /dev/null @@ -1,150 +0,0 @@ -// The MIT License (MIT) -// Copyright (c) 2017 Matt Overby -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. - -#ifndef MCL_BACKTRACKING_H -#define MCL_BACKTRACKING_H - -#include "Problem.hpp" - -namespace mcl { -namespace optlib { - - -// -// Old reliable backtracking Armijo -// -template -class Backtracking { -public: - typedef Eigen::Matrix VecX; - - static inline Scalar search(int verbose, int max_iters, Scalar decrease, const VecX &x, const VecX &p, Problem &problem, Scalar alpha0) { - - // First things first, check descent norm - const Scalar t_eps = std::numeric_limits::epsilon(); - if( p.norm() <= t_eps ){ return decrease; } - - const Scalar tau = 0.7; - Scalar alpha = alpha0; - VecX grad; - if( DIM == Eigen::Dynamic ){ grad = VecX::Zero(x.rows()); } - Scalar fx0 = problem.gradient(x, grad); - Scalar gtp = grad.dot(p); - - int iter = 0; - for( ; iter < max_iters; ++iter ){ - Scalar fxa = problem.value(x + alpha*p); - Scalar fx0_fxa = fx0 + alpha*decrease*gtp; // Armijo condition I - if( fxa <= fx0_fxa ){ break; } // sufficient decrease - alpha *= tau; - } - - if( iter >= max_iters ){ - if( verbose > 0 ){ printf("Backtracking::search Error: Reached max_iters\n"); } - return -1; - } - - return alpha; - } - -}; // end class Backtracking - - -// -// Backtracking Armijo with cubic interpolation -// -template -class BacktrackingCubic { -public: - typedef Eigen::Matrix VecX; - - static inline Scalar search(int verbose, int max_iters, Scalar decrease, const VecX &x, const VecX &p, Problem &problem, Scalar alpha0) { - - // First things first, check descent norm - const Scalar t_eps = std::numeric_limits::epsilon(); - if( p.norm() <= t_eps ){ return decrease; } - - Scalar alpha = alpha0; - VecX grad; - if( DIM == Eigen::Dynamic ){ grad = VecX::Zero(x.rows()); } - Scalar fx0 = problem.gradient(x, grad); - Scalar gtp = grad.dot(p); - Scalar fxp = fx0; - Scalar alphap = alpha; - - int iter = 0; - for( ; iter < max_iters; ++iter ){ - Scalar fxa = problem.value(x + alpha*p); - Scalar fx0_fxa = fx0 + alpha*decrease*gtp; // Armijo condition I - if( fxa <= fx0_fxa ){ break; } // sufficient decrease - - Scalar alpha_tmp = iter == 0 ? - ( gtp / (2.0 * (fx0 + gtp - fxa)) ) : - cubic( fx0, gtp, fxa, alpha, fxp, alphap ); - fxp = fxa; - alphap = alpha; - alpha = range( alpha_tmp, 0.1*alpha, 0.5*alpha ); - } - - if( iter >= max_iters ){ - if( verbose > 0 ){ printf("BacktrackingCubic::search Error: Reached max_iters\n"); } - return -1; - } - - return alpha; - } - -private: - static inline Scalar range( Scalar alpha, Scalar low, Scalar high ){ - if( alpha < low ){ return low; } - else if( alpha > high ){ return high; } - return alpha; - } - - // Cubic interpolation - // fx0 = f(x0) - // gtp = f'(x0)^T p - // fxa = f(x0 + alpha*p) - // alpha = step length - // fxp = previous fxa - // alphap = previous alpha - static inline Scalar cubic( Scalar fx0, Scalar gtp, Scalar fxa, Scalar alpha, Scalar fxp, Scalar alphap ){ - typedef Eigen::Matrix Vec2; - typedef Eigen::Matrix Mat2; - - Scalar mult = 1.0 / ( alpha*alpha * alphap*alphap * (alpha-alphap) ); - Mat2 A; - A(0,0) = alphap*alphap; A(0,1) = -alpha*alpha; - A(1,0) = -alphap*alphap*alphap; A(1,1) = alpha*alpha*alpha; - Vec2 B; - B[0] = fxa - fx0 - alpha*gtp; B[1] = fxp - fx0 - alphap*gtp; - Vec2 r = mult * A * B; - if( std::abs(r[0]) <= 0.0 ){ return -gtp / (2.0*r[1]); } // if quadratic - Scalar d = std::sqrt( r[1]*r[1] - 3.0*r[0]*gtp ); // discrim - return (-r[1] + d) / (3.0*r[0]); - } - -}; // end class BacktrackingCubic - -} // ns optlib -} // ns mcl - -#endif diff --git a/deps/mcloptlib/include/MCL/LBFGS.hpp b/deps/mcloptlib/include/MCL/LBFGS.hpp deleted file mode 100644 index bd06032..0000000 --- a/deps/mcloptlib/include/MCL/LBFGS.hpp +++ /dev/null @@ -1,158 +0,0 @@ -// The MIT License (MIT) -// Copyright (c) 2017 University of Minnesota -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. - -#ifndef MCL_LBFGS_H -#define MCL_LBFGS_H - -#include "Minimizer.hpp" - -namespace mcl { -namespace optlib { - -// L-BFGS implementation based on Nocedal & Wright Numerical Optimization book (Section 7.2) -// DIM = dimension of the problem -// M = history window -// -// Original Author: Ioannis Karamouzas -// -template -class LBFGS : public Minimizer { -private: - typedef Eigen::Matrix VecX; - typedef Eigen::Matrix MatM; - typedef Eigen::Matrix VecM; - -public: - bool show_denom_warning; // Print out warning for zero denominators - - LBFGS() : show_denom_warning(false) { - this->m_settings.max_iters = 50; - show_denom_warning = this->m_settings.verbose > 0 ? true : false; - } - - // Returns number of iterations used - int minimize(Problem &problem, VecX &x){ - - MatM s, y; - VecM alpha, rho; - VecX grad, q, grad_old, x_old, x_last; - - if( DIM==Eigen::Dynamic ){ - int dim = x.rows(); - s = MatM::Zero(dim,M); - y = MatM::Zero(dim,M); - alpha = VecM::Zero(M); - rho = VecM::Zero(M); - grad = VecX::Zero(dim); - q = VecX::Zero(dim); - grad_old = VecX::Zero(dim); - x_old = VecX::Zero(dim); - x_last = VecX::Zero(dim); - } - - problem.gradient(x, grad); - - Scalar gamma_k = 1.0; - Scalar alpha_init = 1.0; - - int global_iter = 0; - int max_iters = this->m_settings.max_iters; - int verbose = this->m_settings.verbose; - - for( int k=0; k= 0; --i){ - rho(i) = 1.0 / ((s.col(i)).dot(y.col(i))); - alpha(i) = rho(i)*(s.col(i)).dot(q); - q = q - alpha(i)*y.col(i); - } - - // L-BFGS second - loop recursion - q = gamma_k*q; - for(int i = 0; i < iter; ++i){ - Scalar beta = rho(i)*q.dot(y.col(i)); - q = q + (alpha(i) - beta)*s.col(i); - } - - // is there a descent - Scalar dir = q.dot(grad); - if(dir <= 0 ){ - q = grad; - max_iters -= k; - k = 0; - alpha_init = std::min(1.0, 1.0 / grad.template lpNorm() ); - } - - Scalar rate = this->linesearch(x, -q, problem, alpha_init); - - if( rate <= 0 ){ - if( verbose > 0 ){ printf("LBFGS::minimize: Failure in linesearch\n"); } - return Minimizer::FAILURE; - } - - x_last = x; - x -= rate * q; - if( problem.converged(x_last,x,grad) ){ break; } - - problem.gradient(x,grad); - VecX s_temp = x - x_old; - VecX y_temp = grad - grad_old; - - // update the history - if(k < M){ - s.col(k) = s_temp; - y.col(k) = y_temp; - } - else { - s.leftCols(M - 1) = s.rightCols(M - 1).eval(); - s.rightCols(1) = s_temp; - y.leftCols(M - 1) = y.rightCols(M - 1).eval(); - y.rightCols(1) = y_temp; - } - - Scalar denom = y_temp.dot(y_temp); - if( std::abs(denom) <= 0 ){ - if( show_denom_warning ){ - printf("LBFGS::minimize Warning: Encountered a zero denominator\n"); - } - break; - } - gamma_k = s_temp.dot(y_temp) / denom; - alpha_init = 1.0; - - } - - return global_iter; - - } // end minimize -}; - -} -} - -#endif diff --git a/deps/mcloptlib/include/MCL/Minimizer.hpp b/deps/mcloptlib/include/MCL/Minimizer.hpp deleted file mode 100644 index 70aa7f8..0000000 --- a/deps/mcloptlib/include/MCL/Minimizer.hpp +++ /dev/null @@ -1,113 +0,0 @@ -// The MIT License (MIT) -// Copyright (c) 2017 Matt Overby -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. - -#ifndef MCL_MINIMIZER_H -#define MCL_MINIMIZER_H - -#include "Problem.hpp" -#include "Backtracking.hpp" -#include "MoreThuente.hpp" -#include "WolfeBisection.hpp" -#include - -namespace mcl { -namespace optlib { - -// The different line search methods currently implemented -enum class LSMethod { - None = 0, // use step length = 1, not recommended - MoreThuente, // TODO test this one for correctness - Backtracking, // basic backtracking with sufficient decrease - BacktrackingCubic, // backtracking with cubic interpolation - WeakWolfeBisection // slow -}; - -// Trust region subproblem method (see TrustRegion.hpp) -enum class TRMethod { - CauchyPoint, - DogLeg -}; - -// -// Base class for optimization algs -// -template -class Minimizer { -public: - typedef Eigen::Matrix VecX; - static const int FAILURE = -1; // returned by minimize if an error is encountered - - struct Settings { - int verbose; // higher = more printouts - int max_iters; // usually changed by derived constructors - int ls_max_iters; // max line search iters - Scalar ls_decrease; // sufficient decrease param - LSMethod ls_method; // see LSMethod (above) - TRMethod tr_method; // see TRMethod (above) - - Settings() : verbose(0), max_iters(100), - ls_max_iters(100000), ls_decrease(1e-4), - ls_method(LSMethod::BacktrackingCubic), - tr_method(TRMethod::DogLeg) - {} - } m_settings; - - // - // Performs optimization - // - virtual int minimize(Problem &problem, VecX &x) = 0; - - -protected: - - // Line search method/options can be changed through m_settings. - Scalar linesearch(const VecX &x, const VecX &p, Problem &prob, double alpha0) const { - double alpha = alpha0; - int mi = m_settings.ls_max_iters; - int v = m_settings.verbose; - Scalar sd = m_settings.ls_decrease; - switch( m_settings.ls_method ){ - default:{ - alpha = Backtracking::search(v, mi, sd, x, p, prob, alpha0); - } break; - case LSMethod::None: { alpha = 1.0; } break; - case LSMethod::MoreThuente: { - alpha = MoreThuente::search(x, p, prob, alpha0); - } break; - case LSMethod::Backtracking: { - alpha = Backtracking::search(v, mi, sd, x, p, prob, alpha0); - } break; - case LSMethod::BacktrackingCubic: { - alpha = BacktrackingCubic::search(v, mi, sd, x, p, prob, alpha0); - } break; - case LSMethod::WeakWolfeBisection: { - alpha = WolfeBisection::search(v, mi, x, p, prob, alpha0); - } break; - } - return alpha; - } // end do linesearch - -}; // class minimizer - -} // ns optlib -} // ns mcl - -#endif diff --git a/deps/mcloptlib/include/MCL/MoreThuente.hpp b/deps/mcloptlib/include/MCL/MoreThuente.hpp deleted file mode 100644 index dd3ff15..0000000 --- a/deps/mcloptlib/include/MCL/MoreThuente.hpp +++ /dev/null @@ -1,330 +0,0 @@ -// The MIT License (MIT) -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. -// -// From https://github.com/PatWie/CppNumericalSolvers -// - -#ifndef MCL_MORETHUENTE_H -#define MCL_MORETHUENTE_H - -#include "Problem.hpp" - -namespace mcl { -namespace optlib { - -template -class MoreThuente { -private: - typedef Eigen::Matrix VectorX; - -public: - - static inline Scalar search(const VectorX &x, const VectorX &p, Problem &problem, Scalar alpha0){ - Scalar alpha = alpha0; - cvsrch(problem, x, alpha, p); - return alpha; - } - - static void cvsrch(Problem &problem, const VectorX &x0, Scalar &stp, const VectorX &s) { - int info = 0; - int infoc = 1; - const Scalar xtol = 1e-15; - const Scalar ftol = 1e-4; - const Scalar gtol = 1e-2; - const Scalar stpmin = 1e-15; - const Scalar stpmax = 1e15; - const Scalar xtrapf = 4; - const int maxfev = 20; - int nfev = 0; - int dim = x0.rows(); - - VectorX g; - if( DIM == Eigen::Dynamic ){ g = VectorX::Zero(dim); } - else{ g.setZero(); } - - Scalar f = problem.gradient(x0, g); - Scalar dginit = g.dot(s); - if (dginit >= 0.0) { - // no descent direction - return; - } - - bool brackt = false; - bool stage1 = true; - - Scalar finit = f; - Scalar dgtest = ftol * dginit; - Scalar width = stpmax - stpmin; - Scalar width1 = 2 * width; - VectorX x = x0.eval(); - - Scalar stx = 0.0; - Scalar fx = finit; - Scalar dgx = dginit; - Scalar sty = 0.0; - Scalar fy = finit; - Scalar dgy = dginit; - - Scalar stmin = 0.0; - Scalar stmax = 0.0; - - const int max_iters = 100000; - int iter = 0; - for( ; iter= stmax))) - || (nfev >= maxfev - 1 ) || (infoc == 0) - || (brackt & (stmax - stmin <= xtol * stmax))) { - stp = stx; - } - - // test new point - x = x0 + stp * s; - f = problem.gradient(x, g); - nfev++; - Scalar dg = g.dot(s); - Scalar ftest1 = finit + stp * dgtest; - - // all possible convergence tests - if ((brackt & ((stp <= stmin) | (stp >= stmax))) | (infoc == 0)) - info = 6; - - if ((stp == stpmax) & (f <= ftest1) & (dg <= dgtest)) - info = 5; - - if ((stp == stpmin) & ((f > ftest1) | (dg >= dgtest))) - info = 4; - - if (nfev >= maxfev) - info = 3; - - if (brackt & (stmax - stmin <= xtol * stmax)) - info = 2; - - if ((f <= ftest1) & (fabs(dg) <= gtol * (-dginit))) - info = 1; - - // terminate when convergence reached - if (info != 0) - return; - - if (stage1 & (f <= ftest1) & (dg >= std::min(ftol, gtol)*dginit)) - stage1 = false; - - if (stage1 & (f <= fx) & (f > ftest1)) { - Scalar fm = f - stp * dgtest; - Scalar fxm = fx - stx * dgtest; - Scalar fym = fy - sty * dgtest; - Scalar dgm = dg - dgtest; - Scalar dgxm = dgx - dgtest; - Scalar dgym = dgy - dgtest; - - cstep( stx, fxm, dgxm, sty, fym, dgym, stp, fm, dgm, brackt, stmin, stmax, infoc); - - fx = fxm + stx * dgtest; - fy = fym + sty * dgtest; - dgx = dgxm + dgtest; - dgy = dgym + dgtest; - } else { - // this is ugly and some variables should be moved to the class scope - cstep( stx, fx, dgx, sty, fy, dgy, stp, f, dg, brackt, stmin, stmax, infoc); - } - - if (brackt) { - if (fabs(sty - stx) >= 0.66 * width1) - stp = stx + 0.5 * (sty - stx); - - width1 = width; - width = fabs(sty - stx); - } - - } // end while true - - if( iter == max_iters ){ - throw std::runtime_error("MoreThuente::linesearch Error: Reached max_iter"); - } - - return; - } - - static void cstep(Scalar& stx, Scalar& fx, Scalar& dx, Scalar& sty, Scalar& fy, Scalar& dy, Scalar& stp, - Scalar& fp, Scalar& dp, bool& brackt, Scalar& stpmin, Scalar& stpmax, int& info) { - info = 0; - bool bound = false; - - // Check the input parameters for errors. - if ((brackt & ((stp <= std::min(stx, sty) ) | (stp >= std::max(stx, sty)))) | (dx * (stp - stx) >= 0.0) - | (stpmax < stpmin)) { - return; - } - - Scalar sgnd = dp * (dx / fabs(dx)); - - Scalar stpf = 0; - Scalar stpc = 0; - Scalar stpq = 0; - - if (fp > fx) { - info = 1; - bound = true; - Scalar theta = 3. * (fx - fp) / (stp - stx) + dx + dp; - Scalar s = std::max(theta, std::max(dx, dp)); - Scalar gamma = s * sqrt((theta / s) * (theta / s) - (dx / s) * (dp / s)); - if (stp < stx) - gamma = -gamma; - - Scalar p = (gamma - dx) + theta; - Scalar q = ((gamma - dx) + gamma) + dp; - Scalar r = p / q; - stpc = stx + r * (stp - stx); - stpq = stx + ((dx / ((fx - fp) / (stp - stx) + dx)) / 2.) * (stp - stx); - if (fabs(stpc - stx) < fabs(stpq - stx)) - stpf = stpc; - else - stpf = stpc + (stpq - stpc) / 2; - - brackt = true; - } else if (sgnd < 0.0) { - info = 2; - bound = false; - Scalar theta = 3 * (fx - fp) / (stp - stx) + dx + dp; - Scalar s = std::max(theta, std::max(dx, dp)); - Scalar gamma = s * sqrt((theta / s) * (theta / s) - (dx / s) * (dp / s)); - if (stp > stx) - gamma = -gamma; - - Scalar p = (gamma - dp) + theta; - Scalar q = ((gamma - dp) + gamma) + dx; - Scalar r = p / q; - stpc = stp + r * (stx - stp); - stpq = stp + (dp / (dp - dx)) * (stx - stp); - if (fabs(stpc - stp) > fabs(stpq - stp)) - stpf = stpc; - else - stpf = stpq; - - brackt = true; - } else if (fabs(dp) < fabs(dx)) { - info = 3; - bound = 1; - Scalar theta = 3 * (fx - fp) / (stp - stx) + dx + dp; - Scalar s = std::max(theta, std::max( dx, dp)); - Scalar gamma = s * sqrt(std::max(static_cast(0.), (theta / s) * (theta / s) - (dx / s) * (dp / s))); - if (stp > stx) - gamma = -gamma; - - Scalar p = (gamma - dp) + theta; - Scalar q = (gamma + (dx - dp)) + gamma; - Scalar r = p / q; - if ((r < 0.0) & (gamma != 0.0)) { - stpc = stp + r * (stx - stp); - } else if (stp > stx) { - stpc = stpmax; - } else { - stpc = stpmin; - } - stpq = stp + (dp / (dp - dx)) * (stx - stp); - if (brackt) { - if (fabs(stp - stpc) < fabs(stp - stpq)) { - stpf = stpc; - } else { - stpf = stpq; - } - } else { - if (fabs(stp - stpc) > fabs(stp - stpq)) { - stpf = stpc; - } else { - stpf = stpq; - } - } - } else { - info = 4; - bound = false; - if (brackt) { - Scalar theta = 3 * (fp - fy) / (sty - stp) + dy + dp; - Scalar s = std::max(theta, std::max(dy, dp)); - Scalar gamma = s * sqrt((theta / s) * (theta / s) - (dy / s) * (dp / s)); - if (stp > sty) - gamma = -gamma; - - Scalar p = (gamma - dp) + theta; - Scalar q = ((gamma - dp) + gamma) + dy; - Scalar r = p / q; - stpc = stp + r * (sty - stp); - stpf = stpc; - } else if (stp > stx) - stpf = stpmax; - else { - stpf = stpmin; - } - } - - if (fp > fx) { - sty = stp; - fy = fp; - dy = dp; - } else { - if (sgnd < 0.0) { - sty = stx; - fy = fx; - dy = dx; - } - stx = stp; - fx = fp; - dx = dp; - } - - stpf = std::min(stpmax, stpf); - stpf = std::max(stpmin, stpf); - stp = stpf; - - if (brackt & bound) { - if (sty > stx) { - stp = std::min(stx + static_cast(0.66) * (sty - stx), stp); - } else { - stp = std::max(stx + static_cast(0.66) * (sty - stx), stp); - } - } - - return; - - } // end cstep - -}; - -} -} - -#endif diff --git a/deps/mcloptlib/include/MCL/Newton.hpp b/deps/mcloptlib/include/MCL/Newton.hpp deleted file mode 100644 index f1632b8..0000000 --- a/deps/mcloptlib/include/MCL/Newton.hpp +++ /dev/null @@ -1,79 +0,0 @@ -// The MIT License (MIT) -// Copyright (c) 2017 Matt Overby -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. - -#ifndef MCL_NEWTON_H -#define MCL_NEWTON_H - -#include "Minimizer.hpp" - -namespace mcl { -namespace optlib { - -template -class Newton : public Minimizer { -private: - typedef Eigen::Matrix VectorX; - typedef Eigen::Matrix MatrixX; - -public: - Newton() { - this->m_settings.max_iters = 20; - } - - int minimize(Problem &problem, VectorX &x){ - - VectorX grad, delta_x, x_last; - if( DIM == Eigen::Dynamic ){ - int dim = x.rows(); - x_last.resize(dim); - grad.resize(dim); - delta_x.resize(dim); - } - - int verbose = this->m_settings.verbose; - int max_iters = this->m_settings.max_iters; - int iter = 0; - for( ; iter < max_iters; ++iter ){ - - problem.gradient(x,grad); - problem.solve_hessian(x,grad,delta_x); - - Scalar rate = this->linesearch(x, delta_x, problem, 1.0); - - if( rate <= 0 ){ - if( verbose > 0 ){ printf("Newton::minimize: Failure in linesearch\n"); } - return Minimizer::FAILURE; - } - - x_last = x; - x += rate * delta_x; - if( problem.converged(x_last,x,grad) ){ break; } - } - - return iter; - } - -}; - -} // ns optlib -} // ns mcl - -#endif diff --git a/deps/mcloptlib/include/MCL/NonLinearCG.hpp b/deps/mcloptlib/include/MCL/NonLinearCG.hpp deleted file mode 100644 index 064fd99..0000000 --- a/deps/mcloptlib/include/MCL/NonLinearCG.hpp +++ /dev/null @@ -1,86 +0,0 @@ -// The MIT License (MIT) -// Copyright (c) 2017 Matt Overby -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. - -#ifndef MCL_NONLINEARCG_H -#define MCL_NONLINEARCG_H - -#include "Minimizer.hpp" - -namespace mcl { -namespace optlib { - -template -class NonLinearCG : public Minimizer { -private: - typedef Eigen::Matrix VectorX; - typedef Eigen::Matrix MatrixX; - -public: - NonLinearCG() { - this->m_settings.max_iters = 100; - } - - int minimize(Problem &problem, VectorX &x){ - - VectorX grad, grad_old, p, x_last; - if( DIM == Eigen::Dynamic ){ - int dim = x.rows(); - x_last.setZero(dim); - grad.resize(dim); - grad_old.resize(dim); - p.resize(dim); - } - - int verbose = this->m_settings.verbose; - int max_iters = this->m_settings.max_iters; - int iter=0; - for( ; iterlinesearch(x, p, problem, 1.0); - - if( rate <= 0 ){ - if( verbose > 0 ){ printf("NonLinearCG::minimize: Failure in linesearch\n"); } - return Minimizer::FAILURE; - } - - x_last = x; - x += rate*p; - grad_old = grad; - - if( problem.converged(x_last,x,grad) ){ break; } - } - return iter; - } // end minimize - -}; - -} // ns optlib -} // ns mcl - -#endif diff --git a/deps/mcloptlib/include/MCL/Problem.hpp b/deps/mcloptlib/include/MCL/Problem.hpp deleted file mode 100644 index 54ebc3c..0000000 --- a/deps/mcloptlib/include/MCL/Problem.hpp +++ /dev/null @@ -1,135 +0,0 @@ -// The MIT License (MIT) -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. - -#ifndef MCL_PROBLEM_H -#define MCL_PROBLEM_H - -#if MCL_DEBUG == 1 -#include -#endif - -#include - -namespace mcl { -namespace optlib { - -template -class Problem { -private: - typedef Eigen::Matrix VecX; - typedef Eigen::Matrix MatX; - -public: - // Returns true if the solver has converged - // x0 is the result of the previous iteration - // x1 is the result at the current iteration - // grad is the gradient at the last iteration - virtual bool converged(const VecX &x0, const VecX &x1, const VecX &grad) = 0; - - // Compute just the value - virtual Scalar value(const VecX &x) = 0; - - // Compute the objective value and the gradient - virtual Scalar gradient(const VecX &x, VecX &grad){ - finiteGradient(x, grad); - return value(x); - } - - // Compute hessian - virtual void hessian(const VecX &x, MatX &hessian){ - finiteHessian(x, hessian); - } - - // Solve dx = H^-1 -g (used by Newton's) - virtual void solve_hessian(const VecX &x, const VecX &grad, VecX &dx){ - MatX hess; - if( DIM == Eigen::Dynamic ){ - int dim = x.rows(); - hess = MatX::Zero(dim,dim); - } - hessian(x,hess); // hessian at x_n - - // Going with with high-accurate, low requirements as default factorization for lin-solve - // Copied from https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html - // Method Requirements Spd (sm) Spd (lg) Accuracy - // partialPivLu() Invertible ++ ++ + - // fullPivLu() None - - - +++ - // householderQr() None ++ ++ + - // colPivHouseholderQr() None ++ - +++ - // fullPivHouseholderQr() None - - - +++ - // llt() PD +++ +++ + - // ldlt() P/N SD +++ + ++ - if( DIM == Eigen::Dynamic || DIM > 4 ){ - dx = hess.fullPivLu().solve(-grad); - } - else { - dx = -hess.inverse()*grad; - } - } - - // Gradient with finite differences - inline void finiteGradient(const VecX &x, VecX &grad){ - const int accuracy = 0; // accuracy can be 0, 1, 2, 3 - const Scalar eps = 2.2204e-6; - const std::vector< std::vector > coeff = - { {1, -1}, {1, -8, 8, -1}, {-1, 9, -45, 45, -9, 1}, {3, -32, 168, -672, 672, -168, 32, -3} }; - const std::vector< std::vector > coeff2 = - { {1, -1}, {-2, -1, 1, 2}, {-3, -2, -1, 1, 2, 3}, {-4, -3, -2, -1, 1, 2, 3, 4} }; - const std::vector dd = {2, 12, 60, 840}; - int dim = x.rows(); - if( grad.rows() != dim ){ grad = VecX::Zero(dim); } - else{ grad.setZero(); } - for(int d = 0; d < dim; ++d){ - for (int s = 0; s < 2*(accuracy+1); ++s){ - VecX xx = x.eval(); - xx[d] += coeff2[accuracy][s]*eps; - grad[d] += coeff[accuracy][s]*value(xx); - } - grad[d] /= (dd[accuracy]* eps); - } - } // end finite grad - - // Hessian with finite differences - inline void finiteHessian(const VecX &x, MatX &hess){ - const Scalar eps = std::numeric_limits::epsilon()*10e7; - int dim = x.rows(); - if( hess.rows() != dim || hess.cols() != dim ){ hess = MatX::Zero(dim,dim); } - for(int i = 0; i < dim; ++i){ - for(int j = 0; j < dim; ++j){ - VecX xx = x; - Scalar f4 = value(xx); - xx[i] += eps; - xx[j] += eps; - Scalar f1 = value(xx); - xx[j] -= eps; - Scalar f2 = value(xx); - xx[j] += eps; - xx[i] -= eps; - Scalar f3 = value(xx); - hess(i, j) = (f1 - f2 - f3 + f4) / (eps * eps); - } - } - } // end finite hess -}; - -} -} - -#endif diff --git a/deps/mcloptlib/include/MCL/TrustRegion.hpp b/deps/mcloptlib/include/MCL/TrustRegion.hpp deleted file mode 100644 index 52bc961..0000000 --- a/deps/mcloptlib/include/MCL/TrustRegion.hpp +++ /dev/null @@ -1,208 +0,0 @@ -// The MIT License (MIT) -// Copyright (c) 2017 Matt Overby -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. - -#ifndef MCL_TRUSTREGION_H -#define MCL_TRUSTREGION_H - -#include "Minimizer.hpp" - -namespace mcl { -namespace optlib { - -template -class TrustRegion : public Minimizer { -private: - typedef Eigen::Matrix VecX; - typedef Eigen::Matrix MatX; - -public: - // - // The trust region method operates on an approximate hessian (B). - // To keep the interface simple, Problem::hessian is used to obtain - // this approximation. It shouldn't be much of a problem, since you - // can still use approximations for Newton's (the only other method - // implemented that requires a hessian evaluation). - // - // TODO Replace B with possibly sparse representation - // - TrustRegion() { - this->m_settings.max_iters = 100; - } - - int minimize(Problem &problem, VecX &x){ - - Scalar delta_k = 2.0; // trust region radius - const Scalar delta_max = 8.0; // max trust region radius - const Scalar eta = 0.125; // min reduction ratio allowed (0m_settings.tr_method; - int max_iters = this->m_settings.max_iters; - int verbose = this->m_settings.verbose; - - VecX grad, dx, x_last; - MatX B; // Approximate hessian - if( DIM == Eigen::Dynamic ){ - int dim = x.rows(); - x_last.resize(dim); // last variable - grad.resize(dim); // grad at k - dx.resize(dim); // descent direction - B.resize(dim,dim); - } - - // Init gradient and hessian - Scalar fxk = problem.gradient(x,grad); // gradient and objective - problem.hessian(x,B); // get hessian (or approximation) - problem.solve_hessian(x,grad,dx); // attempt with newtons - - int iter = 0; - for( ; iter < max_iters; ++iter ){ - - // If it's outside the trust region, pick new descent - if( dx.norm() > delta_k ){ - eval_subproblem(m, delta_k, grad, B, dx); - } - - Scalar fxdx = problem.value(x+dx); - - // Compute reduction ratio - Scalar rho_k = eval_reduction(fxk, fxdx, dx, grad, B); - Scalar dx_norm = dx.norm(); - - // Update trust region radius - if( rho_k < 0.25 ){ delta_k = 0.25*delta_k; } - - // Full step, good approximation - else if( rho_k > 0.75 && std::abs(dx_norm-delta_k) <= 0.0 ){ - delta_k = std::min( 2.0*delta_k, delta_max ); - } - - // Take a step, otherwise need to re-eval sub problem - if( rho_k > eta ){ - - x_last = x; - x = x + dx; - - // Only need to compute gradient and hessian - // if x has actually changed. - fxk = problem.gradient(x,grad); // gradient and objective - if( problem.converged(x_last,x,grad) ){ break; } - - // I think I should improve this as to not call both hessian - // and solve_hessian, which likely causes redundant computation. - problem.hessian(x,B); // get hessian (or approximation) - problem.solve_hessian(x,grad,dx); // attempt with newtons - } - - if( std::isnan(rho_k) ){ - if( verbose ){ printf("\n**TrustRegion Error: NaN reduction"); } - return Minimizer::FAILURE; - } - - } - - return iter; - - } // end minimize - -protected: - - // Assumes coeffs size 3 - static inline Scalar max_roots(Scalar a, Scalar b, Scalar c){ - - Scalar d = (b*b) - (4.0*a*c); - if( d > 0 ){ - Scalar sqrt_d = std::sqrt(d); - Scalar r1 = (-b + sqrt_d) / (2.0*a); - Scalar r2 = (-b - sqrt_d) / (2.0*a); - return std::max(r1,r2); - } - // Should I do something with the real/imaginary parts? - // Scalar real_part = -b/(2.0*a); - // Scalar imaginary = std::sqrt(-d)/(2.0*a); - throw std::runtime_error("TrustRegion Error: Problem in quadratic roots"); - return 0; - } - - // fxk = objective at x_k - // fxdx = objective at x_k + dx - // dx = descent direction - // grad_k = gradient at x_k - // B_k = hessian guess. - static inline Scalar eval_reduction( Scalar fxk, Scalar fxdx, const VecX &dx, - const VecX &grad_k, const MatX &B_k ){ - // rho = ( f(x) - f(x-dx) ) / ( model(0) - model(dx) ) - // with model = f(x) + dx^T grad + 0.5 dx^T B dx - Scalar num = fxk - fxdx; - Scalar denom = fxk - ( fxk + dx.dot(grad_k) + 0.5 * dx.dot( B_k * dx ) ); - return num/denom; - } - - static inline void eval_subproblem( - const TRMethod &m, Scalar delta_k, - const VecX &grad, const MatX &B, - VecX &dx ){ - - Scalar gTBg = grad.dot(B*grad); - - switch( m ){ - - // Generally requires a large number of outer solver iterations - case TRMethod::CauchyPoint: { - - Scalar grad_norm = grad.norm(); - Scalar tau = 1.0; - if( gTBg > 0.0 ){ tau = std::min( 1.0, std::pow(grad_norm,3) / (delta_k*gTBg) ); } - dx = ( -tau * delta_k / grad_norm ) * grad; - - } break; - - // Uses steepest descent if possible - case TRMethod::DogLeg: { - - Scalar gTg = grad.dot(grad); - VecX dx_U = ( -gTg / gTBg ) * grad; - Scalar dx_U_norm = dx_U.norm(); - - // Use steepest descent - if( dx_U_norm >= delta_k ){ dx = delta_k/dx_U_norm * dx_U; } - else{ - // Compute tau and update descent - VecX dx_C = dx - dx_U; - Scalar dx_C_norm = dx_C.norm(); - Scalar tau = max_roots( // Ax^2 + Bx + c - dx_C_norm*dx_C_norm, - 2.0*dx_C.dot(dx_U), - dx_U_norm*dx_U_norm - delta_k*delta_k - ); - dx = dx_U + tau*dx_C; - } - - } break; - - } // end swithc method - - } // end eval sub problem - -}; - -} // ns optlib -} // ns mcl - -#endif diff --git a/deps/mcloptlib/include/MCL/WolfeBisection.hpp b/deps/mcloptlib/include/MCL/WolfeBisection.hpp deleted file mode 100644 index b9e6a06..0000000 --- a/deps/mcloptlib/include/MCL/WolfeBisection.hpp +++ /dev/null @@ -1,110 +0,0 @@ -// The MIT License (MIT) -// Copyright (c) 2018 Matt Overby -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. - -#ifndef MCL_WOLFEBISECTION_H -#define MCL_WOLFEBISECTION_H - -#include "Problem.hpp" - -namespace mcl { -namespace optlib { - -// Bisection method for Weak Wolfe conditions -template -class WolfeBisection { -private: - // Strong wolfe conditions: the armijo rule and a stronger curvature condition - // alpha = step length - // fx_ap = f(x + alpha p) - // fx = f(x) - // pT_gx = p^T ( grad f(x) ) - // pT_gx_ap = p^T ( grad f(x + alpha p) ) - static inline bool strong_wolfe( Scalar alpha, - Scalar fx, Scalar pT_gx, - Scalar fx_ap, Scalar pT_gx_ap, - Scalar wolfe_c1, Scalar wolfe_c2 ){ - if( !(fx_ap <= fx + wolfe_c1 * alpha * pT_gx) ){ return false; } // armijo rule - return std::abs( pT_gx_ap ) <= wolfe_c2 * std::abs( pT_gx ); - } - -public: - typedef Eigen::Matrix VecX; - typedef Eigen::Matrix MatX; - - static inline Scalar search(int verbose, int max_iters, const VecX &x, const VecX &p, Problem &problem, Scalar alpha0) { - - const Scalar t_eps = std::numeric_limits::epsilon(); - const Scalar wolfe_c1 = 0.0001; - const Scalar wolfe_c2 = 0.8; // should be 0.1 for CG! - const int dim = x.rows(); - double alpha = alpha0; - double alpha_min = 1e-8; - double alpha_max = 1; - - VecX grad0, grad_new; - if( DIM == Eigen::Dynamic ){ - grad0 = VecX::Zero(dim); - grad_new = VecX::Zero(dim); - } - Scalar fx0 = problem.gradient(x, grad0); - const Scalar gtp = grad0.dot(p); - bool min_set = false; - - int iter = 0; - for( ; iter < max_iters; ++iter ){ - - // Should we stop iterating? - if( std::abs(alpha_max-alpha_min) <= t_eps ){ break; } - - // Step halfway - alpha = ( alpha_max + alpha_min ) * 0.5; - grad_new.setZero(); - Scalar fx_ap = problem.gradient(x + alpha*p, grad_new); - Scalar gt_ap = grad_new.dot( p ); - - // Check the wolfe conditions - bool happy_wolfe = strong_wolfe( alpha, fx0, gtp, fx_ap, gt_ap, wolfe_c1, wolfe_c2 ); - if( happy_wolfe ){ - alpha_min = alpha; - min_set = true; - } - else { alpha_max = alpha; } - - } // end bs iters - - if( iter == max_iters ){ - if( verbose > 0 ){ printf("WolfeBisection::linesearch Error: Reached max_iters\n"); } - return -1; - } - if( !min_set ){ - if( verbose > 0 ){ printf("WolfeBisection::linesearch Error: LS blocked\n"); } - return -1; - } - - return alpha_min; - - } -}; - -} -} - -#endif diff --git a/deps/mcloptlib/test/TestProblem.hpp b/deps/mcloptlib/test/TestProblem.hpp deleted file mode 100644 index c4c8a4f..0000000 --- a/deps/mcloptlib/test/TestProblem.hpp +++ /dev/null @@ -1,130 +0,0 @@ -// The MIT License (MIT) -// Copyright (c) 2017 Matt Overby -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. - -#include "MCL/Problem.hpp" - -// min |Ax-b| -class DynProblem : public mcl::optlib::Problem { -public: - typedef Eigen::Matrix VectorX; - typedef Eigen::Matrix MatrixX; - - MatrixX A; - VectorX b; - DynProblem( int dim_ ){ - - // Test on random SPD - A = MatrixX::Random(dim_,dim_); - A = A.transpose() * A; - A = A + MatrixX::Identity(dim_,dim_); - - b = VectorX::Random(dim_); - } - - int dim() const { return b.rows(); } - bool converged(const VectorX &x0, const VectorX &x1, const VectorX &grad){ - - // Check sizes of input - int m_dim = dim(); - if( x0.rows() != m_dim ){ - throw std::runtime_error("Error in Problem::converged: x0 wrong dimension"); - } - if( x1.rows() != m_dim ){ - throw std::runtime_error("Error in Problem::converged: x1 wrong dimension"); - } - if( grad.rows() != m_dim ){ - throw std::runtime_error("Error in Problem::converged: gradient wrong dimension"); - } - - return grad.norm() < 1e-10 || (x0-x1).norm() < 1e-10; - } - - double value(const VectorX &x){ - - // Check sizes of input - int m_dim = dim(); - if( x.rows() != m_dim ){ - throw std::runtime_error("Error in Problem::value: x wrong dimension"); - } - - return (A*x-b).norm(); - } - - double gradient(const VectorX &x, VectorX &grad){ - - // Check sizes of input - int m_dim = dim(); - if( x.rows() != m_dim ){ - throw std::runtime_error("Error in Problem::gradient: x wrong dimension"); - } - if( grad.rows() != m_dim ){ - throw std::runtime_error("Error in Problem::gradient: gradient wrong dimension"); - } - - grad = A*x-b; return value(x); - } - void hessian(const VectorX &x, MatrixX &hess){ - - // Check sizes of input - int m_dim = dim(); - if( x.rows() != m_dim ){ - throw std::runtime_error("Error in Problem::hessian: x wrong dimension"); - } - if( hess.rows() != m_dim || hess.cols() != m_dim ){ - throw std::runtime_error("Error in Problem::hessian: hessian wrong dimension"); - } - hess = A; - } - - void solve_hessian(const VectorX &x, const VectorX &grad, VectorX &dx){ - - // Check sizes of input - int m_dim = dim(); - if( x.rows() != m_dim ){ - throw std::runtime_error("Error in Problem::solve_hessian: x wrong dimension"); - } - if( dx.rows() != m_dim ){ - throw std::runtime_error("Error in Problem::solve_hessian: dx wrong dimension"); - } - if( grad.rows() != m_dim ){ - throw std::runtime_error("Error in Problem::solve_hessian: gradient wrong dimension"); - } - - // Check to make sure base class function works as expected - Problem::solve_hessian(x,grad,dx); - } -}; - -class Rosenbrock : public mcl::optlib::Problem { -public: - typedef Eigen::Matrix VectorX; - bool converged(const VectorX &x0, const VectorX &x1, const VectorX &grad){ - (void)(x1); (void)(x0); - return grad.norm() < 1e-10; - } - double value(const VectorX &x){ - double a = 1.0 - x[0]; - double b = x[1] - x[0]*x[0]; - return a*a + b*b*100.0; - } - // Test finite diff as well I guess -}; - diff --git a/deps/mcloptlib/test/testSolvers.cpp b/deps/mcloptlib/test/testSolvers.cpp deleted file mode 100644 index 4f047bf..0000000 --- a/deps/mcloptlib/test/testSolvers.cpp +++ /dev/null @@ -1,210 +0,0 @@ -// The MIT License (MIT) -// Copyright (c) 2017 Matt Overby -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. - -#include -#include "TestProblem.hpp" -#include "MCL/LBFGS.hpp" -#include "MCL/NonLinearCG.hpp" -#include "MCL/Newton.hpp" -#include "MCL/TrustRegion.hpp" -#include - -using namespace mcl::optlib; -typedef std::shared_ptr< Minimizer > MinPtr2; // rb -typedef std::shared_ptr< Minimizer > MinPtrD; // linear - -bool test_linear( std::vector &solvers, std::vector &names ){ - - std::cout << "\nTest linear:" << std::endl; - typedef Eigen::Matrix VecX; - bool success = true; - - // Using multiple dimensions in the linear case, which changes - // how Newton's performs a solve. - std::vector test_dims = { 4, 16, 64 }; - - int n_test_dims = test_dims.size(); - for( int i=0; im_settings.max_iters = 1; } - else { solvers[i]->m_settings.max_iters = 100; } - solvers[i]->m_settings.verbose = 1; - VecX x = VecX::Zero(dim); - solvers[i]->minimize( cp, x ); - for( int i=0; i 1e-3 ){ - std::cerr << "(" << names[i] << ") Failed to minimize: |Ax-b| = " << rn << std::endl; - curr_success = false; - } - - if( curr_success ){ std::cout << "(" << names[i] << ") Linear (" << dim << "): Success" << std::endl; } - else{ success = false; } - - } // end loop solvers - - } // end loop test dims - - return success; - -} - - -bool test_rb( std::vector &solvers, std::vector &names ){ - - std::cout << "\nTest Rosenbrock:" << std::endl; - Rosenbrock rb; // Dim = 2, also tests finite gradient/hessian - bool success = true; - - int n_solvers = solvers.size(); - for( int i=0; im_settings.max_iters = 1000; - solvers[i]->m_settings.verbose = 1; - Eigen::Vector2d x = Eigen::Vector2d::Zero(); - solvers[i]->minimize( rb, x ); - - for( int i=0; i<2; ++i ){ - if( std::isnan(x[i]) || std::isinf(x[i]) ){ - std::cerr << "(" << names[i] << ") Bad values in x: " << x[i] << std::endl; - curr_success = false; - } - } - - double rn = (Eigen::Vector2d(1,1) - x).norm(); - if( rn > 1e-8 ){ - std::cerr << "(" << names[i] << ") Failed to minimize: Rosenbrock = " << rn << std::endl; - curr_success = false; - } - - if( curr_success ){ std::cout << "(" << names[i] << ") Rosenbrock: Success" << std::endl; } - else{ success = false; } - } - - return success; -} - -// Test what happens when the energy is already minimized -bool test_zero( std::vector &solvers, std::vector &names ){ - - std::cout << "\nTest zero energy:" << std::endl; - typedef Eigen::Matrix VecX; - bool success = true; - int dim = 2; - - // High/Low dimensions, linear, and Eigen::Dynamic. - // The solvers should work, since the Vecs/Mats are resized at run time. - DynProblem cp(dim); - - int n_solvers = solvers.size(); - for( int i=0; im_settings.max_iters = 1; - solvers[i]->m_settings.verbose = 1; - VecX x = cp.A.inverse() * cp.b; - solvers[i]->minimize( cp, x ); - - for( int i=0; i 1e-10 ){ - std::cerr << "(" << names[i] << ") Failed to minimize: |Ax-b| = " << rn << std::endl; - curr_success = false; - } - - if( curr_success ){ std::cout << "(" << names[i] << ") Linear (" << dim << "): Success" << std::endl; } - else{ success = false; } - - } // end loop solvers - - return success; -} - - -int main(int argc, char *argv[] ){ - srand(100); - std::vector< std::string > names; - std::vector< MinPtr2 > min2; - std::vector< MinPtrD > minD; - - std::string mode = "all"; - if( argc == 2 ){ mode = std::string(argv[1]); } - - if( mode=="lbfgs" || mode=="all" ){ - min2.emplace_back( std::make_shared< LBFGS >( LBFGS() ) ); - minD.emplace_back( std::make_shared< LBFGS >( LBFGS() ) ); - names.emplace_back( "lbfgs" ); - } - if( mode=="cg" || mode=="all" ){ - min2.emplace_back( std::make_shared< NonLinearCG >( NonLinearCG() ) ); - minD.emplace_back( std::make_shared< NonLinearCG >( NonLinearCG() ) ); - names.emplace_back( "cg" ); - } - if( mode=="newton" || mode=="all" ){ - min2.emplace_back( std::make_shared< Newton >( Newton() ) ); - minD.emplace_back( std::make_shared< Newton >( Newton() ) ); - names.emplace_back( "newton" ); - } - if( mode=="trustregion" || mode=="all" ){ - min2.emplace_back( std::make_shared< TrustRegion >( TrustRegion() ) ); - minD.emplace_back( std::make_shared< TrustRegion >( TrustRegion() ) ); - names.emplace_back( "trustregion" ); - } - - bool success = true; - success &= test_linear( minD, names ); - success &= test_rb( min2, names ); - success &= test_zero( minD, names ); - if( success ){ - std::cout << "\nSUCCESS!" << std::endl; - return EXIT_SUCCESS; - } - else{ std::cout << "\n**FAILURE!" << std::endl; } - return EXIT_FAILURE; -} - - - diff --git a/src/Interp.hpp b/src/Interp.hpp index 4269055..9520d60 100644 --- a/src/Interp.hpp +++ b/src/Interp.hpp @@ -20,7 +20,7 @@ #ifndef INTERP_HPP #define INTERP_HPP 1 -#include +#include #include namespace mpm diff --git a/src/LBFGS.hpp b/src/LBFGS.hpp new file mode 100644 index 0000000..839ebb3 --- /dev/null +++ b/src/LBFGS.hpp @@ -0,0 +1,440 @@ +// Copyright Matt Overby 2021. +// Distributed under the MIT License. +// From: https://github.com/mattoverby/mclgeom + +#ifndef MCL_LBFGS_HPP +#define MCL_LBFGS_HPP 1 + +#include +#include +#include + +namespace mcl +{ + +// L-BFGS implementation based on Nocedal & Wright Numerical Optimization book (Section 7.2) +// adapted from code by Ioannis Karamouzas. +// +// Function pointers are meant to be used with lambdas, e.g. +// LBFGS lbfgs; +// lbfgs.gradient = [&](const MatrixXd &x, MatrixXd &g)->Scalar { return ... }; +// double obj = lbfgs.minimize(x); +// +// If the g arg is not sized, don't compute gradient. This is to avoid +// extra work when only the objective is needed (line search) and to +// avoid redundant code. Example: +// lbfgs.gradient = [&](const MatrixXd &x, MatrixXd &g)->Scalar +// { +// objective = ... +// if (g.rows() == x.rows()) +// { +// g = ... (reuse computation for objective) +// } +// return objective; +// }; +// +template +class LBFGS +{ +public: + typedef typename MatrixType::Scalar Scalar; + + struct Options + { + int min_iters; + int max_iters; + Scalar abs_tol; // absolute tol if converged(...) not set + Scalar rel_tol; // relative tol if converged(...) not set + int M; // history window size + Scalar gamma; // init Hessian = gamma * I + Options() : + min_iters(0), + max_iters(100), + abs_tol(1e-5), + rel_tol(0), + M(6), + gamma(1) + {} + } options; + + LBFGS(); + + // Output from the last call to minimize(...) + int iters() const { return num_iters; } + Scalar gamma() const { return gamma_k; } + + // Resizes buffers and sets to zero. + // Called during minimize(...) ONLY if there is a change in dof or M. + // Otherwise it's assumed you're picking up where you left off + // on the previous call to minimize(...) + void reset(int rows, int cols); + + // Required: + // computes objective value and gradient + // obj = gradient(x, g) + // If the g arg is not sized, don't compute gradient + std::function gradient; + + // Optional: + // Returns true if the solver should exit, default uses ||g|| converged; + + // Optional: + // Linesearch function, default uses bracketing weak wolfe (slow!) + // obj_k1 = (x, grad, descent, alpha) + // Returns new objective value and updates both x AND gradient + std::function linesearch; + + // Optional: + // Filter descent direction, p = B(p) + // Otherwise p = gamma_k * p is used. + std::function filter; + + // Calls initialize(x) once and iterate(x) until converged + Scalar minimize(MatrixType& x); + + // Initialize the solver + void initialize(MatrixType& x); + + // Take an iteration + // Returns objective + Scalar iterate(MatrixType& x); + + // i.e. bisection with weak Wolfe conditions + Scalar bracketing_weakwolfe( + MatrixType& x, + MatrixType& grad, + const MatrixType& p, + Scalar &alpha) const; + + // Used if converged not set + // Returns true if: + // grad.norm <= abs_tol + // or + // grad.norm() <= rel_tol * x.norm() + bool default_converged( + Scalar curr_obj, + const MatrixType& xprev, + const MatrixType& x, + const MatrixType& grad) const; + + Scalar inner( + const MatrixType &a, + const MatrixType &b) const; + +protected: + bool initialized; + int num_iters, max_iters, k; + Scalar gamma_k, obj_0, obj_k; + std::vector s; + std::vector y; + Eigen::Matrix alpha; + Eigen::Matrix rho; + MatrixType grad; + MatrixType q; + MatrixType descent; + MatrixType grad_old; + MatrixType x_old; + MatrixType x_last; + MatrixType s_temp; + MatrixType y_temp; + +}; // end class LBFGS + +// +// Implementation +// + +template +LBFGS::LBFGS() : + initialized(false), + num_iters(0), + max_iters(0), + k(0), + gamma_k(1), + obj_0(0), + obj_k(0) + {} + +template +void LBFGS::reset(int rows, int cols) +{ + using namespace Eigen; + num_iters = 0; + max_iters = 0; + k = 0; + gamma_k = 1; + obj_0 = std::numeric_limits::max(); + obj_k = std::numeric_limits::max(); + int M = options.M; + s = std::vector(M, MatrixType::Zero(rows,cols)); + y = std::vector(M, MatrixType::Zero(rows,cols)); + alpha = VectorXd::Zero(M); + rho = VectorXd::Zero(M); + grad = MatrixType::Zero(rows, cols); + q = MatrixType::Zero(rows, cols); // inv descent + descent = q; + grad_old = MatrixType::Zero(rows, cols); + x_old = MatrixType::Zero(rows, cols); + x_last = MatrixType::Zero(rows, cols); + s_temp = MatrixType::Zero(rows, cols); + y_temp = MatrixType::Zero(rows, cols); +} + +// Returns number of iterations used +template +typename LBFGS::Scalar +LBFGS::minimize(MatrixType &x) +{ + initialize(x); + + // Did we start at the initializer? + if (num_iters >= options.min_iters && + converged(obj_k,x_last,x,grad)) + { + num_iters = 1; + return obj_k; + } + + for (; k +void LBFGS::initialize(MatrixType& x) +{ + initialized = true; + + if (gradient == nullptr) { + throw std::runtime_error("no gradient function"); + } + + if (converged == nullptr) + { + using namespace std::placeholders; + converged = std::bind(&LBFGS::default_converged, this, _1, _2, _3, _4); + } + + if (linesearch == nullptr) + { + using namespace std::placeholders; + linesearch = std::bind(&LBFGS::bracketing_weakwolfe, this, _1, _2, _3, _4); + } + + // Resize/resize variables? + if (alpha.rows() != options.M || + grad.rows() != x.rows() || + grad.cols() != x.cols()) { + reset(x.rows(), x.cols()); + } + + num_iters = 0; + max_iters = std::max(options.min_iters, options.max_iters); + gamma_k = options.gamma; + obj_0 = gradient(x, grad); + obj_k = obj_0; + k = 0; +} + +template +typename LBFGS::Scalar +LBFGS::iterate(MatrixType& x) +{ + if (!initialized) { + throw std::runtime_error("not initialized"); + } + + x_old = x; + grad_old = grad; + q = grad; + num_iters++; + + // + // Two-loop recursion + // + { + // L-BFGS first - loop recursion + int iter = std::min(options.M, k); + for(int i = iter - 1; i >= 0; --i) + { + Scalar denom = inner(s[i], y[i]); + if (std::abs(denom) <= 0.0) + { + rho(i) = 0; + alpha(i) = 0; + continue; + } + rho(i) = 1.0 / denom; + alpha(i) = rho(i)*inner(s[i], q); + q -= alpha(i) * y[i]; + } + + if (filter != nullptr) { filter(q); } + else { q = gamma_k*q; } + + // L-BFGS second - loop recursion + for(int i = 0; i < iter; ++i) + { + Scalar beta = rho(i)*inner(q, y[i]); + q += (alpha(i) - beta)*s[i]; + } + } + + // + // Perform step + // + { + // If our hess approx is bad and we start going in + // the wrong direction, restart memory + Scalar step_size = 1.0; + Scalar dir = inner(q, grad); + if (dir <= 0) + { + q = grad; + max_iters -= k; // Restart memory + k = 0; + step_size = std::min(1.0, 1.0 / grad.template lpNorm() ); + } + + // We've hit local minima, we have to exit + if (q.squaredNorm() <= 0.0) { + return obj_k; + } + + descent = -q; + x_last = x; + obj_k = linesearch(x, grad, descent, step_size); + if (num_iters >= options.min_iters && converged(obj_k,x_last,x,grad)) { + return obj_k; + } + } + + // + // Correction term + // + { + s_temp = x - x_old; + y_temp = grad - grad_old; + + // update the history + if (k < options.M) + { + s[k] = s_temp; + y[k] = y_temp; + } + else + { + for (int i=0; i 0.0) + { + gamma_k = inner(s_temp, y_temp) / denom; + } + } + + return obj_k; +} + +template +typename LBFGS::Scalar +LBFGS::bracketing_weakwolfe( + MatrixType& x, + MatrixType& g, + const MatrixType &drt, + Scalar &step) const +{ + using namespace Eigen; + + Scalar c1 = 10e-4; + Scalar c2 = 0.9; + if (step <= 0.0) { step = 1.0; } + int x_cols = x.cols(); + + g = MatrixType::Zero(x.rows(), x_cols); + MatrixType xp = x; + const Scalar fx_init = gradient(x,g); + Scalar fx = fx_init; + + Scalar dg_init = inner(g, drt); + if (dg_init > 0) { + throw std::runtime_error("direction increases objective"); + } + + const Scalar test_decr = c1 * dg_init; + Scalar lower = 0; + Scalar upper = std::numeric_limits::infinity(); + int maxiter = 200; + int iter = 0; + for (iter = 0; iter < maxiter; iter++) + { + x = xp + step * drt; + fx = gradient(x, g); + if (fx > fx_init + step * test_decr) { // Armijo rule + upper = step; + } + else + { + Scalar dg = inner(g, drt); + if(dg < c2 * dg_init){ // Weak wolfe + lower = step; + } else { + break; // both met + } + } + step = std::isinf(upper) ? 2*step : lower/2 + upper/2; + } + return fx; + +} // end linesearch + +template +bool LBFGS::default_converged( + Scalar curr_obj, + const MatrixType& xprev, + const MatrixType& x, + const MatrixType& g) const +{ + (void)(curr_obj); + (void)(xprev); + Scalar gnorm = g.norm(); + if (gnorm <= options.abs_tol) { + return true; + } + Scalar xnorm = x.norm(); + if (gnorm < options.rel_tol * xnorm) { + return true; + } + return false; +} + +template +typename LBFGS::Scalar +LBFGS::inner( + const MatrixType &a, + const MatrixType &b) const +{ + int cols = std::min(a.cols(), b.cols()); + Scalar dot = 0; + for (int i=0; i namespace mpm { diff --git a/src/Particle.hpp b/src/Particle.hpp index 1dd9693..6be8521 100644 --- a/src/Particle.hpp +++ b/src/Particle.hpp @@ -22,6 +22,7 @@ #include #include +#include #include "Interp.hpp" namespace mpm diff --git a/src/Solver.cpp b/src/Solver.cpp index 847fb31..13be7ad 100644 --- a/src/Solver.cpp +++ b/src/Solver.cpp @@ -22,30 +22,6 @@ namespace mpm { -// The system as represented by an optimization problem for generic solvers -class Objective : public mcl::optlib::Problem -{ -public: - Objective(Solver *solver_) : solver(solver_) {} - Solver *solver; - - double value(const Eigen::VectorXd &v) - { - Eigen::VectorXd grad; - return gradient(v,grad); - } - - bool converged(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1, const Eigen::VectorXd &grad) - { - if (grad.norm() < 1e-2) { return true; } - if ((x0-x1).norm() < 1e-2) { return true; } - return false; - } - - double gradient(const Eigen::VectorXd &v, Eigen::VectorXd &grad); - -}; // end class Objective - bool Solver::initialize() { // Fill a sphere with particles @@ -78,8 +54,6 @@ bool Solver::initialize() std::cout << "Num Particles: " << m_particles.size() << std::endl; std::cout << "Num Nodes: " << m_grid.size() << std::endl; - // Set up the solver - optimizer.m_settings.ls_method = mcl::optlib::LSMethod::Backtracking; return true; } @@ -169,7 +143,7 @@ void Solver::explicit_solve() } // end compute forces -double Objective::gradient(const Eigen::VectorXd &v, Eigen::VectorXd &grad) +double Solver::gradient(const Eigen::VectorXd &v, Eigen::VectorXd &grad) { using namespace Eigen; double tot_energy = 0.0; @@ -178,8 +152,8 @@ double Objective::gradient(const Eigen::VectorXd &v, Eigen::VectorXd &grad) // Compute Particle Deformation Gradients for new grid velocities // #pragma omp parallel for reduction (+:tot_energy) - for(int i=0; im_particles.size(); ++i){ - Particle *p = solver->m_particles[i]; + for(int i=0; iget_deform_grad(); Matrix3d vel_grad; vel_grad.fill(0.f); @@ -190,7 +164,7 @@ double Objective::gradient(const Eigen::VectorXd &v, Eigen::VectorXd &grad) MPM::make_grid_loop(p, grid_nodes, wip, dwip); for(int j=0; jm_grid[ grid_nodes[j] ]->active_idx; + int idx = m_grid[ grid_nodes[j] ]->active_idx; Eigen::Vector3d curr_v(v[idx*3+0], v[idx*3+1], v[idx*3+2]); vel_grad += (MPM::timestep_s*curr_v) * dwip[j].transpose(); } @@ -209,9 +183,9 @@ double Objective::gradient(const Eigen::VectorXd &v, Eigen::VectorXd &grad) // Compute energy gradient // #pragma omp parallel for reduction (+:tot_energy) - for(int i=0; iactive_grid.size(); ++i) + for(int i=0; iactive_grid[i]; + GridNode *node = active_grid[i]; Eigen::Vector3d curr_v(v[i*3+0], v[i*3+1], v[i*3+2]); Eigen::Vector3d momentum_grad = node->m * (curr_v - node->v); @@ -249,8 +223,13 @@ void Solver::implicit_solve() } // Minimize - Objective obj(this); - optimizer.minimize(obj, v); + optimizer.gradient = [&](const VectorXd &x, VectorXd &g)->double + { + double objective = gradient(x, g); + return objective; + }; + + optimizer.minimize(v); // Copy solver results back to grid #pragma omp parallel for diff --git a/src/Solver.hpp b/src/Solver.hpp index 127cf4c..53ebdbc 100644 --- a/src/Solver.hpp +++ b/src/Solver.hpp @@ -21,7 +21,7 @@ #define MPM_SOLVER_HPP 1 #include "MPM.hpp" -#include "MCL/LBFGS.hpp" +#include "LBFGS.hpp" namespace mpm { @@ -41,7 +41,7 @@ class Solver std::vector active_grid; // resized each time step std::vector m_particles; void get_vertices(Eigen::MatrixXd &X) const; - mcl::optlib::LBFGS optimizer; + mcl::LBFGS optimizer; // // Settings @@ -61,6 +61,9 @@ class Solver // Returns true on success. bool step(float screen_dt); + // Gradient of the objective. + double gradient(const Eigen::VectorXd &v, Eigen::VectorXd &grad); + private: double compute_timestep(double screen_dt); @@ -72,7 +75,7 @@ class Solver void explicit_solve(); void implicit_solve(); -}; // end class system +}; // end class solver } // end namespace mpm diff --git a/test/solver.cpp b/test/solver.cpp new file mode 100644 index 0000000..ddffc6f --- /dev/null +++ b/test/solver.cpp @@ -0,0 +1,9 @@ +#include + +// placeholder +int main(int argc, char *argv[]) +{ + (void)(argc); + (void)(argv); + return EXIT_SUCCESS; +} \ No newline at end of file From fa862dbb319e6cabaf5c4878f1a3fb1c80de2189 Mon Sep 17 00:00:00 2001 From: mattoverby Date: Wed, 8 Oct 2025 14:19:27 -0500 Subject: [PATCH 2/2] mv sphere to test --- .github/workflows/build_and_test.yml | 30 ++++++++++++++++ CMakeLists.txt | 51 +++++++--------------------- test/solver.cpp | 3 ++ {src => test}/sphere.cpp | 0 4 files changed, 45 insertions(+), 39 deletions(-) create mode 100644 .github/workflows/build_and_test.yml rename {src => test}/sphere.cpp (100%) diff --git a/.github/workflows/build_and_test.yml b/.github/workflows/build_and_test.yml new file mode 100644 index 0000000..2d619d7 --- /dev/null +++ b/.github/workflows/build_and_test.yml @@ -0,0 +1,30 @@ +name: Build and Test + +on: + push: + branches: [ master ] + pull_request: + branches: [ master ] + +jobs: + build_and_test: + runs-on: ubuntu-latest + + steps: + - name: clone + uses: actions/checkout@v4 # Action to check out your repository + + - name: dependencies + run: | + sudo apt-get update + sudo apt-get install -y build-essential cmake + + - name: build + working-directory: ${{ github.workspace }} + run: | + cmake -DMPM_TESTING_ONLY=ON . + make -j + + - name: test + working-directory: ${{ github.workspace }} + run: make test \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 8fc6554..179767f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,42 +1,18 @@ -# Copyright (c) 2016 Matt Overby -# -# MPM-OPTIMIZATION Uses the BSD 2-Clause License (http://www.opensource.org/licenses/BSD-2-Clause) -# Redistribution and use in source and binary forms, with or without modification, are -# permitted provided that the following conditions are met: -# 1. Redistributions of source code must retain the above copyright notice, this list of -# conditions and the following disclaimer. -# 2. Redistributions in binary form must reproduce the above copyright notice, this list -# of conditions and the following disclaimer in the documentation and/or other materials -# provided with the distribution. -# THIS SOFTWARE IS PROVIDED "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -# ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF MINNESOTA, DULUTH OR CONTRIBUTORS BE -# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, -# OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER -# IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY -# OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +cmake_minimum_required(VERSION 3.27) +set(CMAKE_EXPORT_COMPILE_COMMANDS ON) -cmake_minimum_required(VERSION 3.0) project(mpm C CXX) +set(CMAKE_CXX_STANDARD 17) +list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake") -set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH}) -set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") -set(CMAKE_BUILD_TYPE Release) - -add_definitions(-DMPM_SRC_DIR="${CMAKE_CURRENT_SOURCE_DIR}") +# Options option(MPM_TESTING_ONLY "Only compile lib and tests" OFF) -############################################################ -# -# Libraries -# -############################################################ - # Libigl also includes Eigen, etc include(libigl) # OpenMP +# (TO-DO: replace with TBB) find_package(OpenMP) if (OPENMP_FOUND) set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") @@ -44,6 +20,7 @@ if (OPENMP_FOUND) add_definitions(-DOMP_NESTED) endif() +# Simulation library set(MPM_SRC ${CMAKE_CURRENT_SOURCE_DIR}/src/LBFGS.hpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Solver.hpp @@ -52,24 +29,20 @@ set(MPM_SRC ${CMAKE_CURRENT_SOURCE_DIR}/src/MPM.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Interp.hpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Particle.hpp) - - add_library(mpm_optimization ${MPM_SRC}) target_link_libraries(mpm_optimization PUBLIC igl::core) -############################################################ -# -# Binaries -# -############################################################ - +# Interactive viewer if (NOT MPM_TESTING_ONLY) igl_include(glfw) - add_executable(sphere ${CMAKE_CURRENT_SOURCE_DIR}/src/sphere.cpp) + add_executable(sphere ${CMAKE_CURRENT_SOURCE_DIR}/test/sphere.cpp) target_link_libraries(sphere PUBLIC mpm_optimization igl::glfw) + target_include_directories(sphere PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/src) endif() +# Unit testing (TO-DO) enable_testing() add_executable(test_solver ${CMAKE_CURRENT_SOURCE_DIR}/test/solver.cpp) target_link_libraries(test_solver PUBLIC mpm_optimization) +target_include_directories(test_solver PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/src) add_test(NAME TestSolver COMMAND test_solver) \ No newline at end of file diff --git a/test/solver.cpp b/test/solver.cpp index ddffc6f..aa06e89 100644 --- a/test/solver.cpp +++ b/test/solver.cpp @@ -1,8 +1,11 @@ #include +#include "Solver.hpp" // placeholder int main(int argc, char *argv[]) { + mpm::Solver solver; + (void)(argc); (void)(argv); return EXIT_SUCCESS; diff --git a/src/sphere.cpp b/test/sphere.cpp similarity index 100% rename from src/sphere.cpp rename to test/sphere.cpp