From 1c7117c6a0a0a2b66d959e442e41429888c62461 Mon Sep 17 00:00:00 2001 From: Iason Date: Mon, 23 Nov 2020 11:06:45 +0200 Subject: [PATCH] Initial commit --- CMakeLists.txt | 74 +++ cmake/DownloadProject.CMakeLists.cmake.in | 17 + cmake/DownloadProject.cmake | 182 ++++++ cmake/FindLIBIGL.cmake | 34 ++ src/gradientDescent.h | 709 ++++++++++++++++++++++ src/main.cpp | 135 ++++ src/reducedmodeloptimizer.cpp | 415 +++++++++++++ src/reducedmodeloptimizer.hpp | 57 ++ 8 files changed, 1623 insertions(+) create mode 100644 CMakeLists.txt create mode 100644 cmake/DownloadProject.CMakeLists.cmake.in create mode 100644 cmake/DownloadProject.cmake create mode 100644 cmake/FindLIBIGL.cmake create mode 100644 src/gradientDescent.h create mode 100644 src/main.cpp create mode 100644 src/reducedmodeloptimizer.cpp create mode 100644 src/reducedmodeloptimizer.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..fec2faa --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,74 @@ +cmake_minimum_required(VERSION 2.8) +project(ReducedModelOptimization) + +#Add the project cmake scripts to the module path +list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake) + +#Download external dependencies NOTE: If the user has one of these libs it shouldn't be downloaded again. +include(${CMAKE_MODULE_PATH}/DownloadProject.cmake) +if (CMAKE_VERSION VERSION_LESS 3.2) + set(UPDATE_DISCONNECTED_IF_AVAILABLE "") +else() + set(UPDATE_DISCONNECTED_IF_AVAILABLE "UPDATE_DISCONNECTED 1") +endif() + +##vcglib devel branch +download_project(PROJ vcglib_devel + GIT_REPOSITORY https://github.com/cnr-isti-vclab/vcglib.git + GIT_TAG devel + ${UPDATE_DISCONNECTED_IF_AVAILABLE} +) + +##matplot++ lib +download_project(PROJ MATPLOTPLUSPLUS + GIT_REPOSITORY https://github.com/alandefreitas/matplotplusplus + GIT_TAG master + ${UPDATE_DISCONNECTED_IF_AVAILABLE} + ) +add_subdirectory(${MATPLOTPLUSPLUS_SOURCE_DIR}) + +#Polyscope +download_project(PROJ POLYSCOPE + GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git + GIT_TAG master + ${UPDATE_DISCONNECTED_IF_AVAILABLE} +) +add_subdirectory(${POLYSCOPE_SOURCE_DIR}) + +#bobyqa-cpp +download_project(PROJ BOBYQA + GIT_REPOSITORY https://github.com/elsid/bobyqa-cpp.git + GIT_TAG master + ${UPDATE_DISCONNECTED_IF_AVAILABLE} +) +add_subdirectory(${BOBYQA_SOURCE_DIR}) +message(STATUS "BOBYQA bin dir:${BOBYQA_BINARY_DIR}") + +##Eigen 3 NOTE: Eigen is required on the system the code is ran +find_package(Eigen3 3.3 REQUIRED) + +#Add the project sources +file(GLOB SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/src/*.hpp ${CMAKE_CURRENT_SOURCE_DIR}/src/*.cpp) + +#OpenMP +find_package(OpenMP) + +set(MYSOURCESDIR "/home/iason/Coding/Libraries/MySources") +file(GLOB MYSOURCES ${MYSOURCESDIR}/*.hpp ${MYSOURCESDIR}/*.cpp) + +add_executable(${PROJECT_NAME} ${SOURCES} ${MYSOURCES} ) + +target_include_directories(${PROJECT_NAME} + PRIVATE ${vcglib_devel_SOURCE_DIR} + PRIVATE ${MYSOURCESDIR} + PRIVATE ${BOBYQA_SOURCE_DIR}/include/ + ) + +#Use C++17 +if(MSVC) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /std:c++20") +else(MSVC) + target_compile_features(${PROJECT_NAME} PUBLIC cxx_std_20) +endif(MSVC) + +target_link_libraries(${PROJECT_NAME} polyscope Eigen3::Eigen OpenMP::OpenMP_CXX matplot bobyqa_shared) diff --git a/cmake/DownloadProject.CMakeLists.cmake.in b/cmake/DownloadProject.CMakeLists.cmake.in new file mode 100644 index 0000000..89be4fd --- /dev/null +++ b/cmake/DownloadProject.CMakeLists.cmake.in @@ -0,0 +1,17 @@ +# Distributed under the OSI-approved MIT License. See accompanying +# file LICENSE or https://github.com/Crascit/DownloadProject for details. + +cmake_minimum_required(VERSION 2.8.2) + +project(${DL_ARGS_PROJ}-download NONE) + +include(ExternalProject) +ExternalProject_Add(${DL_ARGS_PROJ}-download + ${DL_ARGS_UNPARSED_ARGUMENTS} + SOURCE_DIR "${DL_ARGS_SOURCE_DIR}" + BINARY_DIR "${DL_ARGS_BINARY_DIR}" + CONFIGURE_COMMAND "" + BUILD_COMMAND "" + INSTALL_COMMAND "" + TEST_COMMAND "" +) diff --git a/cmake/DownloadProject.cmake b/cmake/DownloadProject.cmake new file mode 100644 index 0000000..e300f42 --- /dev/null +++ b/cmake/DownloadProject.cmake @@ -0,0 +1,182 @@ +# Distributed under the OSI-approved MIT License. See accompanying +# file LICENSE or https://github.com/Crascit/DownloadProject for details. +# +# MODULE: DownloadProject +# +# PROVIDES: +# download_project( PROJ projectName +# [PREFIX prefixDir] +# [DOWNLOAD_DIR downloadDir] +# [SOURCE_DIR srcDir] +# [BINARY_DIR binDir] +# [QUIET] +# ... +# ) +# +# Provides the ability to download and unpack a tarball, zip file, git repository, +# etc. at configure time (i.e. when the cmake command is run). How the downloaded +# and unpacked contents are used is up to the caller, but the motivating case is +# to download source code which can then be included directly in the build with +# add_subdirectory() after the call to download_project(). Source and build +# directories are set up with this in mind. +# +# The PROJ argument is required. The projectName value will be used to construct +# the following variables upon exit (obviously replace projectName with its actual +# value): +# +# projectName_SOURCE_DIR +# projectName_BINARY_DIR +# +# The SOURCE_DIR and BINARY_DIR arguments are optional and would not typically +# need to be provided. They can be specified if you want the downloaded source +# and build directories to be located in a specific place. The contents of +# projectName_SOURCE_DIR and projectName_BINARY_DIR will be populated with the +# locations used whether you provide SOURCE_DIR/BINARY_DIR or not. +# +# The DOWNLOAD_DIR argument does not normally need to be set. It controls the +# location of the temporary CMake build used to perform the download. +# +# The PREFIX argument can be provided to change the base location of the default +# values of DOWNLOAD_DIR, SOURCE_DIR and BINARY_DIR. If all of those three arguments +# are provided, then PREFIX will have no effect. The default value for PREFIX is +# CMAKE_BINARY_DIR. +# +# The QUIET option can be given if you do not want to show the output associated +# with downloading the specified project. +# +# In addition to the above, any other options are passed through unmodified to +# ExternalProject_Add() to perform the actual download, patch and update steps. +# The following ExternalProject_Add() options are explicitly prohibited (they +# are reserved for use by the download_project() command): +# +# CONFIGURE_COMMAND +# BUILD_COMMAND +# INSTALL_COMMAND +# TEST_COMMAND +# +# Only those ExternalProject_Add() arguments which relate to downloading, patching +# and updating of the project sources are intended to be used. Also note that at +# least one set of download-related arguments are required. +# +# If using CMake 3.2 or later, the UPDATE_DISCONNECTED option can be used to +# prevent a check at the remote end for changes every time CMake is run +# after the first successful download. See the documentation of the ExternalProject +# module for more information. It is likely you will want to use this option if it +# is available to you. Note, however, that the ExternalProject implementation contains +# bugs which result in incorrect handling of the UPDATE_DISCONNECTED option when +# using the URL download method or when specifying a SOURCE_DIR with no download +# method. Fixes for these have been created, the last of which is scheduled for +# inclusion in CMake 3.8.0. Details can be found here: +# +# https://gitlab.kitware.com/cmake/cmake/commit/bdca68388bd57f8302d3c1d83d691034b7ffa70c +# https://gitlab.kitware.com/cmake/cmake/issues/16428 +# +# If you experience build errors related to the update step, consider avoiding +# the use of UPDATE_DISCONNECTED. +# +# EXAMPLE USAGE: +# +# include(DownloadProject) +# download_project(PROJ googletest +# GIT_REPOSITORY https://github.com/google/googletest.git +# GIT_TAG master +# UPDATE_DISCONNECTED 1 +# QUIET +# ) +# +# add_subdirectory(${googletest_SOURCE_DIR} ${googletest_BINARY_DIR}) +# +#======================================================================================== + + +set(_DownloadProjectDir "${CMAKE_CURRENT_LIST_DIR}") + +include(CMakeParseArguments) + +function(download_project) + + set(options QUIET) + set(oneValueArgs + PROJ + PREFIX + DOWNLOAD_DIR + SOURCE_DIR + BINARY_DIR + # Prevent the following from being passed through + CONFIGURE_COMMAND + BUILD_COMMAND + INSTALL_COMMAND + TEST_COMMAND + ) + set(multiValueArgs "") + + cmake_parse_arguments(DL_ARGS "${options}" "${oneValueArgs}" "${multiValueArgs}" ${ARGN}) + + # Hide output if requested + if (DL_ARGS_QUIET) + set(OUTPUT_QUIET "OUTPUT_QUIET") + else() + unset(OUTPUT_QUIET) + message(STATUS "Downloading/updating ${DL_ARGS_PROJ}") + endif() + + # Set up where we will put our temporary CMakeLists.txt file and also + # the base point below which the default source and binary dirs will be. + # The prefix must always be an absolute path. + if (NOT DL_ARGS_PREFIX) + set(DL_ARGS_PREFIX "${CMAKE_BINARY_DIR}") + else() + get_filename_component(DL_ARGS_PREFIX "${DL_ARGS_PREFIX}" ABSOLUTE + BASE_DIR "${CMAKE_CURRENT_BINARY_DIR}") + endif() + if (NOT DL_ARGS_DOWNLOAD_DIR) + set(DL_ARGS_DOWNLOAD_DIR "${DL_ARGS_PREFIX}/${DL_ARGS_PROJ}-download") + endif() + + # Ensure the caller can know where to find the source and build directories + if (NOT DL_ARGS_SOURCE_DIR) + set(DL_ARGS_SOURCE_DIR "${DL_ARGS_PREFIX}/${DL_ARGS_PROJ}-src") + endif() + if (NOT DL_ARGS_BINARY_DIR) + set(DL_ARGS_BINARY_DIR "${DL_ARGS_PREFIX}/${DL_ARGS_PROJ}-build") + endif() + set(${DL_ARGS_PROJ}_SOURCE_DIR "${DL_ARGS_SOURCE_DIR}" PARENT_SCOPE) + set(${DL_ARGS_PROJ}_BINARY_DIR "${DL_ARGS_BINARY_DIR}" PARENT_SCOPE) + + # The way that CLion manages multiple configurations, it causes a copy of + # the CMakeCache.txt to be copied across due to it not expecting there to + # be a project within a project. This causes the hard-coded paths in the + # cache to be copied and builds to fail. To mitigate this, we simply + # remove the cache if it exists before we configure the new project. It + # is safe to do so because it will be re-generated. Since this is only + # executed at the configure step, it should not cause additional builds or + # downloads. + file(REMOVE "${DL_ARGS_DOWNLOAD_DIR}/CMakeCache.txt") + + # Create and build a separate CMake project to carry out the download. + # If we've already previously done these steps, they will not cause + # anything to be updated, so extra rebuilds of the project won't occur. + # Make sure to pass through CMAKE_MAKE_PROGRAM in case the main project + # has this set to something not findable on the PATH. + configure_file("${_DownloadProjectDir}/DownloadProject.CMakeLists.cmake.in" + "${DL_ARGS_DOWNLOAD_DIR}/CMakeLists.txt") + execute_process(COMMAND ${CMAKE_COMMAND} -G "${CMAKE_GENERATOR}" + -D "CMAKE_MAKE_PROGRAM:FILE=${CMAKE_MAKE_PROGRAM}" + . + RESULT_VARIABLE result + ${OUTPUT_QUIET} + WORKING_DIRECTORY "${DL_ARGS_DOWNLOAD_DIR}" + ) + if(result) + message(FATAL_ERROR "CMake step for ${DL_ARGS_PROJ} failed: ${result}") + endif() + execute_process(COMMAND ${CMAKE_COMMAND} --build . + RESULT_VARIABLE result + ${OUTPUT_QUIET} + WORKING_DIRECTORY "${DL_ARGS_DOWNLOAD_DIR}" + ) + if(result) + message(FATAL_ERROR "Build step for ${DL_ARGS_PROJ} failed: ${result}") + endif() + +endfunction() diff --git a/cmake/FindLIBIGL.cmake b/cmake/FindLIBIGL.cmake new file mode 100644 index 0000000..1017008 --- /dev/null +++ b/cmake/FindLIBIGL.cmake @@ -0,0 +1,34 @@ +# - Try to find the LIBIGL library +# Once done this will define +# +# LIBIGL_FOUND - system has LIBIGL +# LIBIGL_INCLUDE_DIR - **the** LIBIGL include directory +if(LIBIGL_FOUND) + return() +endif() + +find_path(LIBIGL_INCLUDE_DIR igl/readOBJ.h + HINTS + ${LIBIGL_DIR} + ENV LIBIGL_DIR + PATHS + ${CMAKE_SOURCE_DIR}/../.. + ${CMAKE_SOURCE_DIR}/.. + ${CMAKE_SOURCE_DIR} + ${CMAKE_SOURCE_DIR}/libigl + ${CMAKE_SOURCE_DIR}/../libigl + ${CMAKE_SOURCE_DIR}/../../libigl + /usr + /usr/local + /usr/local/igl/libigl + PATH_SUFFIXES include +) + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(LIBIGL + "\nlibigl not found --- You can download it using:\n\tgit clone https://github.com/libigl/libigl.git ${CMAKE_SOURCE_DIR}/../libigl" + LIBIGL_INCLUDE_DIR) +mark_as_advanced(LIBIGL_INCLUDE_DIR) + +list(APPEND CMAKE_MODULE_PATH "${LIBIGL_INCLUDE_DIR}/../cmake") +include(libigl) diff --git a/src/gradientDescent.h b/src/gradientDescent.h new file mode 100644 index 0000000..0471eca --- /dev/null +++ b/src/gradientDescent.h @@ -0,0 +1,709 @@ +/* gdcpp.h + * + * Author: Fabian Meyer + * Created On: 12 Jul 2019 + * License: MIT + */ + +#ifndef GDCPP_GDCPP_H_ +#define GDCPP_GDCPP_H_ + +#include +#include +#include +#include +#include + +namespace gdc { +typedef long int Index; + +/** Functor to compute forward differences. + * Computes the gradient of the objective f(x) as follows: + * + * grad(x) = (f(x + eps) - f(x)) / eps + * + * The computation requires len(x) evaluations of the objective. + */ +template class ForwardDifferences { +public: + typedef Eigen::Matrix Vector; + typedef std::function Objective; + +private: + Scalar eps_; + Index threads_; + Objective objective_; + +public: + ForwardDifferences() + : ForwardDifferences(std::sqrt(std::numeric_limits::epsilon())) {} + + ForwardDifferences(const Scalar eps) : eps_(eps), threads_(1), objective_() {} + + void setNumericalEpsilon(const Scalar eps) { eps_ = eps; } + + void setThreads(const Index threads) { threads_ = threads; } + + void setObjective(const Objective &objective) { objective_ = objective; } + + void operator()(const Vector &xval, const Scalar fval, Vector &gradient) { + assert(objective_); + + gradient.resize(xval.size()); +#pragma omp parallel for num_threads(threads_) + for (Index i = 0; i < xval.size(); ++i) { + Vector xvalN = xval; + xvalN(i) += eps_; + Scalar fvalN = objective_(xvalN); + + gradient(i) = (fvalN - fval) / eps_; + } + } +}; + +/** Functor to compute backward differences. + * Computes the gradient of the objective f(x) as follows: + * + * grad(x) = (f(x) - f(x - eps)) / eps + * + * The computation requires len(x) evaluations of the objective. + */ +template class BackwardDifferences { +public: + typedef Eigen::Matrix Vector; + typedef std::function Objective; + +private: + Scalar eps_; + Index threads_; + Objective objective_; + +public: + BackwardDifferences() + : BackwardDifferences(std::sqrt(std::numeric_limits::epsilon())) { + } + + BackwardDifferences(const Scalar eps) + : eps_(eps), threads_(1), objective_() {} + + void setNumericalEpsilon(const Scalar eps) { eps_ = eps; } + + void setThreads(const Index threads) { threads_ = threads; } + + void setObjective(const Objective &objective) { objective_ = objective; } + + void operator()(const Vector &xval, const Scalar fval, Vector &gradient) { + assert(objective_); + + gradient.resize(xval.size()); +#pragma omp parallel for num_threads(threads_) + for (Index i = 0; i < xval.size(); ++i) { + Vector xvalN = xval; + xvalN(i) -= eps_; + Scalar fvalN = objective_(xvalN); + gradient(i) = (fval - fvalN) / eps_; + } + } +}; + +/** Functor to compute central differences. + * Computes the gradient of the objective f(x) as follows: + * + * grad(x) = (f(x + 0.5 eps) - f(x - 0.5 eps)) / eps + * + * The computation requires 2 * len(x) evaluations of the objective. + */ +template struct CentralDifferences { +public: + typedef Eigen::Matrix Vector; + typedef std::function Objective; + +private: + Scalar eps_; + Index threads_; + Objective objective_; + +public: + CentralDifferences() + : CentralDifferences(std::sqrt(std::numeric_limits::epsilon())) {} + + CentralDifferences(const Scalar eps) : eps_(eps), threads_(1), objective_() {} + + void setNumericalEpsilon(const Scalar eps) { eps_ = eps; } + + void setThreads(const Index threads) { threads_ = threads; } + + void setObjective(const Objective &objective) { objective_ = objective; } + + void operator()(const Vector &xval, const Scalar, Vector &gradient) { + assert(objective_); + + Vector fvals(xval.size() * 2); +#pragma omp parallel for num_threads(threads_) + for (Index i = 0; i < fvals.size(); ++i) { + Index idx = i / 2; + Vector xvalN = xval; + if (i % 2 == 0) + xvalN(idx) += eps_ / 2; + else + xvalN(idx) -= eps_ / 2; + + fvals(i) = objective_(xvalN); + } + + gradient.resize(xval.size()); + for (Index i = 0; i < xval.size(); ++i) + gradient(i) = (fvals(i * 2) - fvals(i * 2 + 1)) / eps_; + } +}; + +/** Dummy callback functor, which does nothing. */ +template struct NoCallback { + typedef Eigen::Matrix Vector; + + bool operator()(const Index, const Vector &, const Scalar, + const Vector &) const { + return true; + } +}; + +/** Step size functor, which returns a constant step size. */ +template class ConstantStepSize { +public: + typedef Eigen::Matrix Vector; + typedef std::function Objective; + typedef std::function + FiniteDifferences; + +private: + Scalar stepSize_; + +public: + ConstantStepSize() : ConstantStepSize(0.000000000000001) {} + + ConstantStepSize(const Scalar stepSize) : stepSize_(stepSize) {} + + /** Set the step size returned by this functor. + * @param stepSize step size returned by functor */ + void setStepSize(const Scalar stepSize) { stepSize_ = stepSize; } + + void setObjective(const Objective &) {} + + void setFiniteDifferences(const FiniteDifferences &) {} + + Scalar operator()(const Vector &, const Scalar, const Vector &) { + return stepSize_; + } +}; + +/** Step size functor to compute Barzilai-Borwein (BB) steps. + * The functor can either compute the direct or inverse BB step. + * The steps are computed as follows: + * + * s_k = x_k - x_k-1 k >= 1 + * y_k = grad_k - grad_k-1 k >= 1 + * Direct: stepSize = (s_k^T * s_k) / (y_k^T * s_k) + * Inverse: stepSize = (y_k^T * s_k) / (y_k^T * y_k) + * + * The very first step is computed as a constant. */ +template class BarzilaiBorwein { +public: + typedef Eigen::Matrix Vector; + typedef std::function Objective; + typedef std::function + FiniteDifferences; + + enum class Method { Direct, Inverse }; + +private: + Vector lastXval_; + Vector lastGradient_; + Method method_; + Scalar constStep_; + + Scalar constantStep() const { return constStep_; } + + Scalar directStep(const Vector &xval, const Vector &gradient) { + auto sk = xval - lastXval_; + auto yk = gradient - lastGradient_; + Scalar num = sk.dot(sk); + Scalar denom = sk.dot(yk); + + if (denom == 0) + return 1; + else + return std::abs(num / denom); + } + + Scalar inverseStep(const Vector &xval, const Vector &gradient) { + auto sk = xval - lastXval_; + auto yk = gradient - lastGradient_; + Scalar num = sk.dot(yk); + Scalar denom = yk.dot(yk); + + if (denom == 0) + return 1; + else + return std::abs(num / denom); + } + +public: + BarzilaiBorwein() : BarzilaiBorwein(Method::Inverse, 1) {} + + BarzilaiBorwein(const Method method, const Scalar constStep) + : lastXval_(), lastGradient_(), method_(method), constStep_(constStep) {} + + void setObjective(const Objective &) {} + + void setFiniteDifferences(const FiniteDifferences &) {} + + void setMethod(const Method method) { method_ = method; } + + void setConstStepSize(const Scalar stepSize) { constStep_ = stepSize; } + + Scalar operator()(const Vector &xval, const Scalar, const Vector &gradient) { + Scalar stepSize = 0; + if (lastXval_.size() == 0) { + stepSize = constStep_; + } else { + switch (method_) { + case Method::Direct: + stepSize = directStep(xval, gradient); + break; + case Method::Inverse: + stepSize = inverseStep(xval, gradient); + break; + default: + assert(false); + break; + } + } + + lastGradient_ = gradient; + lastXval_ = xval; + + return stepSize; + } +}; + +/** Step size functor to perform Armijo Linesearch with backtracking. + * The functor iteratively decreases the step size until the following + * conditions are met: + * + * Armijo: f(x - stepSize * grad(x)) <= f(x) - cArmijo * stepSize * grad(x)^T * + * grad(x) + * + * If either condition does not hold the step size is decreased: + * + * stepSize = decrease * stepSize */ +template class ArmijoBacktracking { +public: + typedef Eigen::Matrix Vector; + typedef std::function Objective; + typedef std::function + FiniteDifferences; + +protected: + Scalar decrease_; + Scalar cArmijo_; + Scalar minStep_; + Scalar maxStep_; + Index maxIt_; + Objective objective_; + FiniteDifferences finiteDifferences_; + + Scalar evaluateObjective(const Vector &xval, Vector &gradient) { + gradient.resize(0); + Scalar fval = objective_(xval, gradient); + if (gradient.size() == 0) + finiteDifferences_(xval, fval, gradient); + return fval; + } + + virtual bool computeSecondCondition(const Scalar, const Scalar, const Scalar, + const Vector &, const Vector &) { + return true; + } + +public: + ArmijoBacktracking() + : ArmijoBacktracking(0.8, 1e-4, 1e-20, 1, 0) {} // NOTE: maxStep was 1 + + ArmijoBacktracking(const Scalar decrease, const Scalar cArmijo, + const Scalar minStep, const Scalar maxStep, + const Index iterations) + : decrease_(decrease), cArmijo_(cArmijo), minStep_(minStep), + maxStep_(maxStep), maxIt_(iterations), objective_() { + assert(decrease > 0); + assert(decrease < 1); + assert(cArmijo > 0); + assert(cArmijo < 0.5); + assert(minStep < maxStep); + } + + /** Set the decreasing factor for backtracking. + * Assure that decrease in (0, 1). + * @param decrease decreasing factor */ + void setBacktrackingDecrease(const Scalar decrease) { + assert(decrease > 0); + assert(decrease < 1); + decrease_ = decrease; + } + + /** Set the relaxation constant for the Armijo condition (see class + * description). + * Assure cArmijo in (0, 0.5). + * @param cArmijo armijo constant */ + void setArmijoConstant(const Scalar cArmijo) { + assert(cArmijo > 0); + assert(cArmijo < 0.5); + cArmijo_ = cArmijo; + } + + /** Set the bounds for the step size during linesearch. + * The final step size is guaranteed to be in [minStep, maxStep]. + * @param minStep minimum step size + * @param maxStep maximum step size */ + void setStepBounds(const Scalar minStep, const Scalar maxStep) { + assert(minStep < maxStep); + minStep_ = minStep; + maxStep_ = maxStep; + } + + /** Set the maximum number of iterations. + * Set to 0 or negative for infinite iterations. + * @param iterations maximum number of iterations */ + void setMaxIterations(const Index iterations) { maxIt_ = iterations; } + + void setObjective(const Objective &objective) { objective_ = objective; } + + void setFiniteDifferences(const FiniteDifferences &finiteDifferences) { + finiteDifferences_ = finiteDifferences; + } + + Scalar operator()(const Vector &xval, const Scalar fval, + const Vector &gradient) { + assert(objective_); + assert(finiteDifferences_); + + Scalar stepSize = maxStep_ / decrease_; + Vector gradientN; + Vector xvalN; + Scalar fvalN; + bool armijoCondition = false; + bool secondCondition = false; + + Index iterations = 0; + while ((maxIt_ <= 0 || iterations < maxIt_) && + stepSize * decrease_ >= minStep_ && + !(armijoCondition && secondCondition)) { + stepSize = decrease_ * stepSize; + xvalN = xval - stepSize * gradient; + fvalN = evaluateObjective(xvalN, gradientN); + + armijoCondition = + fvalN <= fval - cArmijo_ * stepSize * gradient.dot(gradient); + secondCondition = + computeSecondCondition(stepSize, fval, fvalN, gradient, gradientN); + + ++iterations; + } + + return stepSize; + } +}; + +/** Step size functor to perform Wolfe Linesearch with backtracking. + * The functor iteratively decreases the step size until the following + * conditions are met: + * + * Armijo: f(x - stepSize * grad(x)) <= f(x) - cArmijo * stepSize * grad(x)^T * + * grad(x) Wolfe: grad(x)^T grad(x - stepSize * grad(x)) <= cWolfe * grad(x)^T * + * grad(x) + * + * If either condition does not hold the step size is decreased: + * + * stepSize = decrease * stepSize */ +template +class WolfeBacktracking : public ArmijoBacktracking { +public: + typedef Eigen::Matrix Vector; + typedef std::function Objective; + typedef std::function + FiniteDifferences; + +protected: + Scalar cWolfe_; + + virtual bool computeSecondCondition(const Scalar, const Scalar, const Scalar, + const Vector &gradient, + const Vector &gradientN) { + return gradient.dot(gradientN) <= cWolfe_ * gradient.dot(gradient); + } + +public: + WolfeBacktracking() : WolfeBacktracking(0.8, 1e-4, 0.9, 1e-20, 1, 0) {} + + WolfeBacktracking(const Scalar decrease, const Scalar cArmijo, + const Scalar cWolfe, const Scalar minStep, + const Scalar maxStep, const Index iterations) + : ArmijoBacktracking(decrease, cArmijo, minStep, maxStep, + iterations), + cWolfe_(cWolfe) { + assert(cWolfe < 1); + assert(cArmijo < cWolfe); + } + + /** Set the wolfe constants for Armijo and Wolfe condition (see class + * description). + * Assure that c1 < c2 < 1 and c1 in (0, 0.5). + * @param c1 armijo constant + * @param c2 wolfe constant */ + void setWolfeConstant(const Scalar cWolfe) { + assert(cWolfe < 1); + cWolfe_ = cWolfe; + } +}; + +/** Step size functor which searches for a step that reduces the function + * value. + * The functor iteratively decreases the step size until the following + * condition is met: + * + * f(x - stepSize * grad) < f(x) + * + * If this condition does not hold the step size is decreased: + * + * stepSize = decrease * stepSize + * + * This functor does not require to compute any gradients and does not use + * finite differences. */ +template class DecreaseBacktracking { +public: + typedef Eigen::Matrix Vector; + typedef std::function Objective; + typedef std::function + FiniteDifferences; + +private: + Scalar decrease_; + Scalar minStep_; + Scalar maxStep_; + Index maxIt_; + Objective objective_; + +public: + DecreaseBacktracking() : DecreaseBacktracking(0.8, 1e-12, 1, 0) {} + + DecreaseBacktracking(const Scalar decrease, const Scalar minStep, + const Scalar maxStep, const Index iterations) + : decrease_(decrease), minStep_(minStep), maxStep_(maxStep), + maxIt_(iterations), objective_() {} + + /** Set the decreasing factor for backtracking. + * Assure that decrease in (0, 1). + * @param decrease decreasing factor */ + void setBacktrackingDecrease(const Scalar decrease) { decrease_ = decrease; } + + /** Set the bounds for the step size during linesearch. + * The final step size is guaranteed to be in [minStep, maxStep]. + * @param minStep minimum step size + * @param maxStep maximum step size */ + void setStepBounds(const Scalar minStep, const Scalar maxStep) { + assert(minStep < maxStep); + minStep_ = minStep; + maxStep_ = maxStep; + } + + /** Set the maximum number of iterations. + * Set to 0 or negative for infinite iterations. + * @param iterations maximum number of iterations */ + void setMaxIterations(const Index iterations) { maxIt_ = iterations; } + + void setObjective(const Objective &objective) { objective_ = objective; } + + void setFiniteDifferences(const FiniteDifferences &) {} + + Scalar operator()(const Vector &xval, const Scalar fval, + const Vector &gradient) { + assert(objective_); + + Scalar stepSize = maxStep_ / decrease_; + Vector xvalN; + Vector gradientN; + Scalar fvalN; + bool improvement = false; + + Index iterations = 0; + while ((maxIt_ <= 0 || iterations < maxIt_) && + stepSize * decrease_ >= minStep_ && !improvement) { + stepSize = decrease_ * stepSize; + xvalN = xval - stepSize * gradient; + fvalN = objective_(xvalN, gradientN); + + improvement = fvalN < fval; + + ++iterations; + } + + return stepSize; + } +}; + +template , + typename Callback = NoCallback, + typename FiniteDifferences = CentralDifferences> +class GradientDescent { +public: + typedef Eigen::Matrix Vector; + + struct Result { + Index iterations; + bool converged; + Scalar fval; + Vector xval; + }; + +protected: + Index maxIt_; + Scalar minGradientLen_; + Scalar minStepLen_; + Scalar momentum_; + Index verbosity_; + Objective objective_; + StepSize stepSize_; + Callback callback_; + FiniteDifferences finiteDifferences_; + + Scalar evaluateObjective(const Vector &xval, Vector &gradient) { + gradient.resize(0); + Scalar fval = objective_(xval, gradient); + if (gradient.size() == 0) + finiteDifferences_(xval, fval, gradient); + return fval; + } + + std::string vector2str(const Vector &vec) const { + std::stringstream ss1; + ss1 << std::fixed << std::showpoint << std::setprecision(16); + std::stringstream ss2; + ss2 << '['; + for (Index i = 0; i < vec.size(); ++i) { + ss1 << vec(i); + ss2 << std::setfill(' ') << std::setw(10) << ss1.str(); + if (i != vec.size() - 1) + ss2 << ' '; + ss1.str(""); + } + ss2 << ']'; + + return ss2.str(); + } + +public: + GradientDescent() + : maxIt_(0), minGradientLen_(static_cast(1e-2)), + minStepLen_(static_cast(1e-6)), momentum_(0), verbosity_(0), + objective_(), stepSize_(), callback_(), finiteDifferences_() {} + + ~GradientDescent() {} + + void setThreads(const Index threads) { + finiteDifferences_.setThreads(threads); + } + + void setNumericalEpsilon(const Scalar eps) { + finiteDifferences_.setNumericalEpsilon(eps); + } + + void setMaxIterations(const Index iterations) { maxIt_ = iterations; } + + void setObjective(const Objective &objective) { objective_ = objective; } + + void setCallback(const Callback &callback) { callback_ = callback; } + + void setMinGradientLength(const Scalar gradientLen) { + minGradientLen_ = gradientLen; + } + + void setMinStepLength(const Scalar stepLen) { minStepLen_ = stepLen; } + + void setStepSize(const StepSize stepSize) { stepSize_ = stepSize; } + + void setMomentum(const Scalar momentum) { momentum_ = momentum; } + + void setVerbosity(const Index verbosity) { verbosity_ = verbosity; } + + Result minimize(const Vector &initialGuess) { + finiteDifferences_.setObjective([this](const Vector &xval) { + Vector tmp; + return this->objective_(xval, tmp); + }); + stepSize_.setObjective([this](const Vector &xval, Vector &gradient) { + return this->objective_(xval, gradient); + }); + stepSize_.setFiniteDifferences( + [this](const Vector &xval, const Scalar fval, Vector &gradient) { + this->finiteDifferences_(xval, fval, gradient); + }); + + Vector xval = initialGuess; + Vector gradient; + Scalar fval; + Scalar gradientLen = minGradientLen_ + 1; + Scalar stepSize; + Vector step = Vector::Zero(xval.size()); + Scalar stepLen = minStepLen_ + 1; + bool callbackResult = true; + + Index iterations = 0; + while ((maxIt_ <= 0 || iterations < maxIt_) && + gradientLen >= minGradientLen_ && stepLen >= minStepLen_ && + callbackResult) { + xval -= step; + fval = evaluateObjective(xval, gradient); + gradientLen = gradient.norm(); + // update step according to step size and momentum + stepSize = stepSize_(xval, fval, gradient); + step = momentum_ * step + (1 - momentum_) * stepSize * gradient; + stepLen = step.norm(); + // evaluate callback an save its result + callbackResult = callback_(iterations, xval, fval, gradient); + + if (verbosity_ > 0) { + std::stringstream ss; + ss << "it=" << std::setfill('0') << std::setw(4) << iterations + << std::fixed << std::showpoint << std::setprecision(20) + << " gradlen=" << gradientLen << " stepsize=" << stepSize + << " steplen=" << stepLen; + + if (verbosity_ > 2) + ss << " callback=" << (callbackResult ? "true" : "false"); + + ss << " fval=" << fval; + + if (verbosity_ > 1) + ss << " xval=" << vector2str(xval); + if (verbosity_ > 2) + ss << " gradient=" << vector2str(gradient); + if (verbosity_ > 3) + ss << " step=" << vector2str(step); + std::cout << ss.str() << std::endl; + } + + ++iterations; + } + + Result result; + result.xval = xval; + result.fval = fval; + result.iterations = iterations; + result.converged = gradientLen < minGradientLen_ || stepLen < minStepLen_; + + return result; + } +}; +} // namespace gdc + +#endif diff --git a/src/main.cpp b/src/main.cpp new file mode 100644 index 0000000..1d03182 --- /dev/null +++ b/src/main.cpp @@ -0,0 +1,135 @@ +#include "beamformfinder.hpp" +#include "edgemesh.hpp" +#include "flatpattern.hpp" +#include "polyscope/curve_network.h" +#include "polyscope/point_cloud.h" +#include "polyscope/polyscope.h" +#include "reducedmodeloptimizer.hpp" +#include "simulationhistoryplotter.hpp" +#include +#include +#include +#include +#include + +int main(int argc, char *argv[]) { + FlatPattern pattern("/home/iason/Models/valid_6777.ply"); + // FlatPattern pattern("/home/iason/Models/simple_beam_paper_example.ply"); + pattern.createFan(); + // pattern.savePly("fannedValid.ply"); + FlatPattern reducedModel("/home/iason/Models/reducedModel.ply"); + reducedModel.createFan(); + + // reducedModel.savePly("fannedReduced.ply"); + + std::unordered_map reducedToFullViMap{ + {1, 3}, {2, 7}, {3, 11}, {4, 15}, {5, 19}, {6, 23}}; + ReducedModelOptimizer optimizer(pattern, reducedModel, reducedToFullViMap, + {}); // 6 + // fan reduced model + // {{1, 3}, {2, 15}}, {}); // 2 fan reduced model + Eigen::VectorXd optimalParameters = optimizer.optimize(); + + // Full model simulation + std::unordered_map> fixedVertices; + // fixedVertices[0] = std::unordered_set{0, 1, 2, 3, 4, 5}; + fixedVertices[3] = std::unordered_set{0, 1, 2, 3, 4, 5}; + // fixedVertices[7] = std::unordered_set{0, 1, 2}; + std::unordered_map nodalForces{ + {15, {0, 0, 2000, 0, 0, 0}}}; + SimulationJob fullModelSimulationJob{ + std::make_shared(pattern), fixedVertices, nodalForces, {}}; + Simulator formFinder; + SimulationResults fullModelResults = + formFinder.executeSimulation(fullModelSimulationJob); + fullModelResults.simulationLabel = "Full Model"; + fullModelResults.draw(fullModelSimulationJob); + // fullModelSimulationJob.draw(); + // double stiffnessFactor = 1; + + // while (true) { + // Reduced model simulation + // SimulationJob reducedModelSimulationJob = + // optimizer.getReducedSimulationJob(fullModelSimulationJob); + // const std::vector stiffnessVector{ + // fullModelSimulationJob.mesh->elements[0].properties.A, + // stiffnessFactor * + // fullModelSimulationJob.mesh->elements[0].properties.J, stiffnessFactor + // * fullModelSimulationJob.mesh->elements[0].properties.I2, + // stiffnessFactor * + // fullModelSimulationJob.mesh->elements[0].properties.I3}; + // for (EdgeIndex ei = 0; ei < reducedModelSimulationJob.mesh->EN(); ei++) { + // BeamFormFinder::Element &e = + // reducedModelSimulationJob.mesh->elements[ei]; e.properties.A = + // 0.00035185827018667374; + // // stifnessVector[0]; + // e.properties.J = // stiffnessVector[1]; + // 7.709325104874406e-08; + // e.properties.I2 = -0.0000015661453308127776; + // // stiffnessVector[2]; + // e.properties.I3 = 3.7099813776947167e-07; // stiffnessVector[3]; + // e.axialConstFactor = e.properties.E * e.properties.A / e.initialLength; + // e.torsionConstFactor = e.properties.G * e.properties.J / + // e.initialLength; e.firstBendingConstFactor = + // 2 * e.properties.E * e.properties.I2 / e.initialLength; + // e.secondBendingConstFactor = + // 2 * e.properties.E * e.properties.I3 / e.initialLength; + // } + + // SimulationResults reducedModelResults = + // formFinder.executeSimulation(reducedModelSimulationJob, true); + // reducedModelResults.simulationLabel = + // "Reduced Model_" + std::to_string(stiffnessFactor); + // reducedModelResults.draw(reducedModelSimulationJob); + // SimulationResultsReporter resultsReporter; + // resultsReporter.reportResults( + // {reducedModelResults}, + // std::filesystem::current_path().append("Results"), + // "_" + std::to_string(stiffnessFactor)); + // if (reducedModelResults.history.numberOfSteps == + // BeamFormFinder::Simulator::maxDRMIterations) { + // break; + // } + // stiffnessFactor *= 1.5; + // } + + // Beam + // VCGEdgeMesh beamMesh; + // beamMesh.loadFromPly("/home/iason/Models/simple_beam_model_2elem_1m.ply"); + // const VertexIndex reducedModelOpposite_vi = + // std::ceil(reducedModelSimulationJob.mesh->VN() / 2.0); + // auto v0 = + // reducedModelSimulationJob.mesh->vert[reducedModelOpposite_vi].cP() - + // reducedModelSimulationJob.mesh->vert[1].cP(); + // auto v1 = beamMesh.vert[2].cP() - beamMesh.vert[0].cP(); + // vcg::Matrix44d R; + // const double rotationTheta = + // std::asin((v0 ^ v1).Norm() / (v0.Norm() * v1.Norm())); + // R.SetRotateRad(-rotationTheta, v0 ^ v1); + // vcg::tri::UpdatePosition::Matrix(beamMesh, R); + // vcg::tri::UpdatePosition::Scale(beamMesh, v0.Norm() / + // v1.Norm()); vcg::tri::UpdatePosition::Translate( + // beamMesh, reducedModelSimulationJob.mesh->vert[1].cP()); + // std::unordered_map> + // beamFixedVertices; + // beamFixedVertices[0] = std::unordered_set{0, 1, 2, 3, 4, 5}; + // std::unordered_map beamNodalForces{ + // {2, Eigen::Vector3d(0, 0, 1000)}}; + // SimulationJob beamSimulationJob{std::make_shared(beamMesh), + // beamFixedVertices, + // beamNodalForces, + // {}}; + // // for (EdgeIndex ei = 0; ei < beamSimulationJob.mesh->EN(); ei++) { + // // BeamFormFinder::Element &e = beamSimulationJob.mesh->elements[ei]; + // // e.properties.A = stifnessVector[0]; + // // e.properties.J = stifnessVector[1]; + // // e.properties.I2 = stifnessVector[2]; + // // e.properties.I3 = stifnessVector[3]; + // // } + // // beamSimulationJob.draw(); + // SimulationResults beamSimulationResults = + // formFinder.executeSimulation(beamSimulationJob); + // beamSimulationResults.simulationLabel = "Beam"; + // beamSimulationResults.draw(beamSimulationJob); + return 0; +} diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp new file mode 100644 index 0000000..8f9fc4d --- /dev/null +++ b/src/reducedmodeloptimizer.cpp @@ -0,0 +1,415 @@ +#include "reducedmodeloptimizer.hpp" +#include "bobyqa.h" +#include "gradientDescent.h" +#include "matplot/matplot.h" +#include "simulationhistoryplotter.hpp" + +const bool draw = true; +size_t numberOfOptimizationRounds{0}; +Simulator simulator; +Eigen::MatrixX3d optimalReducedModelDisplacements; +SimulationJob reducedModelSimulationJob; +std::unordered_map + reducedToFullVertexIndexMap; +matplot::line_handle plotHandle; +std::vector fValueHistory; +Eigen::Vector4d initialParameters; +struct OptimizationCallback { + double operator()(const size_t &iterations, const Eigen::VectorXd &x, + const double &fval, Eigen::VectorXd &gradient) const { + // run simulation + // SimulationResults reducedModelResults = + // simulator.executeSimulation(reducedModelSimulationJob); + // reducedModelResults.draw(reducedModelSimulationJob); + fValueHistory.push_back(fval); + auto xPlot = + matplot::linspace(0, fValueHistory.size(), fValueHistory.size()); + plotHandle = matplot::scatter(xPlot, fValueHistory); + // const std::string plotImageFilename = "objectivePlot.png"; + // matplot::save(plotImageFilename); + // if (numberOfOptimizationRounds % 30 == 0) { + // std::filesystem::copy_file( + // std::filesystem::path(plotImageFilename), + // std::filesystem::path("objectivePlot_copy.png")); + // } + // std::stringstream ss; + // ss << x; + // reducedModelResults.simulationLabel = ss.str(); + // SimulationResultsReporter resultsReporter; + // resultsReporter.reportResults( + // {reducedModelResults}, + // std::filesystem::current_path().append("Results")); + + return true; + } +}; + +struct Objective { + double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const { + assert(x.rows() == 4); + + // drawSimulationJob(simulationJob); + // Set mesh from x + std::shared_ptr reducedModel = + reducedModelSimulationJob.mesh; + for (EdgeIndex ei = 0; ei < reducedModel->EN(); ei++) { + BeamFormFinder::Element &e = reducedModel->elements[ei]; + e.properties.A = initialParameters(0) * x(0); + e.properties.J = initialParameters(1) * x(1); + e.properties.I2 = initialParameters(2) * x(2); + e.properties.I3 = initialParameters(3) * x(3); + e.axialConstFactor = e.properties.E * e.properties.A / e.initialLength; + e.torsionConstFactor = e.properties.G * e.properties.J / e.initialLength; + e.firstBendingConstFactor = + 2 * e.properties.E * e.properties.I2 / e.initialLength; + e.secondBendingConstFactor = + 2 * e.properties.E * e.properties.I3 / e.initialLength; + } + // run simulation + SimulationResults reducedModelResults = + simulator.executeSimulation(reducedModelSimulationJob); + // std::stringstream ss; + // ss << x; + // reducedModelResults.simulationLabel = ss.str(); + // SimulationResultsReporter resultsReporter; + // resultsReporter.reportResults( + // {reducedModelResults}, + // std::filesystem::current_path().append("Results")); + // compute error and return it + double error = 0; + for (const auto reducedFullViPair : reducedToFullVertexIndexMap) { + VertexIndex reducedModelVi = reducedFullViPair.first; + Eigen::Vector3d vertexDisplacement( + reducedModelResults.displacements[reducedModelVi][0], + reducedModelResults.displacements[reducedModelVi][1], + reducedModelResults.displacements[reducedModelVi][2]); + Eigen::Vector3d errorVector = + Eigen::Vector3d( + optimalReducedModelDisplacements.row(reducedModelVi)) - + vertexDisplacement; + error += errorVector.norm(); + } + + return error; + } +}; + +double objective(long n, const double *x) { + Eigen::VectorXd eigenX(n, 1); + for (size_t xi = 0; xi < n; xi++) { + eigenX(xi) = x[xi]; + } + // Eigen::VectorXd emptyGradient; + // return operator()(eigenX, emptyGradient); + + // drawSimulationJob(simulationJob); + // Set mesh from x + std::shared_ptr reducedModel = reducedModelSimulationJob.mesh; + for (EdgeIndex ei = 0; ei < reducedModel->EN(); ei++) { + BeamFormFinder::Element &e = reducedModel->elements[ei]; + e.properties.A = initialParameters(0) * eigenX(0); + e.properties.J = initialParameters(1) * eigenX(1); + e.properties.I2 = initialParameters(2) * eigenX(2); + e.properties.I3 = initialParameters(3) * eigenX(3); + e.axialConstFactor = e.properties.E * e.properties.A / e.initialLength; + e.torsionConstFactor = e.properties.G * e.properties.J / e.initialLength; + e.firstBendingConstFactor = + 2 * e.properties.E * e.properties.I2 / e.initialLength; + e.secondBendingConstFactor = + 2 * e.properties.E * e.properties.I3 / e.initialLength; + } + // run simulation + SimulationResults reducedModelResults = + simulator.executeSimulation(reducedModelSimulationJob); + // std::stringstream ss; + // ss << x; + // reducedModelResults.simulationLabel = ss.str(); + // SimulationResultsReporter resultsReporter; + // resultsReporter.reportResults( + // {reducedModelResults}, + // std::filesystem::current_path().append("Results")); + // compute error and return it + double error = 0; + for (const auto reducedFullViPair : reducedToFullVertexIndexMap) { + VertexIndex reducedModelVi = reducedFullViPair.first; + Eigen::Vector3d vertexDisplacement( + reducedModelResults.displacements[reducedModelVi][0], + reducedModelResults.displacements[reducedModelVi][1], + reducedModelResults.displacements[reducedModelVi][2]); + Eigen::Vector3d errorVector = + Eigen::Vector3d(optimalReducedModelDisplacements.row(reducedModelVi)) - + vertexDisplacement; + error += errorVector.norm(); + } + + return error; +} + +ReducedModelOptimizer::ReducedModelOptimizer( + ConstVCGEdgeMesh &fullModel, ConstVCGEdgeMesh &reducedModel, + const std::unordered_map + &viMap, + const std::vector> &oppositeVertices) + : reducedToFullViMap(viMap) { + reducedToFullVertexIndexMap = viMap; + + // construct fullToReducedViMap + for (const std::pair &viPair : + viMap) { + fullToReducedViMap[viPair.second] = viPair.first; + } + + // construct opposite vertex map + for (const std::pair &viPair : oppositeVertices) { + oppositeVertexMap[viPair.first] = viPair.second; + oppositeVertexMap[viPair.second] = viPair.first; + } + + reducedModelElementalMesh = std::make_shared(reducedModel); + pFullModelElementalMesh = std::make_shared(fullModel); + optimalReducedModelDisplacements.resize(reducedModelElementalMesh->VN(), 3); +} + +void ReducedModelOptimizer::computeReducedModelSimulationJob( + const SimulationJob &simulationJobOfFullModel, + SimulationJob &simulationJobOfReducedModel) { + std::unordered_map> + reducedModelFixedVertices; + for (auto fullModelFixedVertex : simulationJobOfFullModel.fixedVertices) { + reducedModelFixedVertices[fullToReducedViMap.at( + fullModelFixedVertex.first)] = fullModelFixedVertex.second; + } + + std::unordered_map reducedModelNodalForces; + for (auto fullModelNodalForce : + simulationJobOfFullModel.nodalExternalForces) { + reducedModelNodalForces[fullToReducedViMap.at(fullModelNodalForce.first)] = + fullModelNodalForce.second; + } + + simulationJobOfReducedModel = SimulationJob{reducedModelElementalMesh, + reducedModelFixedVertices, + reducedModelNodalForces, + {}}; +} + +SimulationJob ReducedModelOptimizer::getReducedSimulationJob( + const SimulationJob &fullModelSimulationJob) { + SimulationJob reducedModelSimulationJob; + computeReducedModelSimulationJob(fullModelSimulationJob, + reducedModelSimulationJob); + return reducedModelSimulationJob; +} + +void ReducedModelOptimizer::computeDesiredReducedModelDisplacements( + const SimulationResults &fullModelResults, + Eigen::MatrixX3d &optimaldisplacementsOfReducedModel) { + for (auto reducedFullViPair : reducedToFullVertexIndexMap) { + const VertexIndex fullModelVi = reducedFullViPair.second; + Vector6d fullModelViDisplacements = + fullModelResults.displacements[fullModelVi]; + optimaldisplacementsOfReducedModel.row(reducedFullViPair.first) = + Eigen::Vector3d(fullModelViDisplacements[0], + fullModelViDisplacements[1], + fullModelViDisplacements[2]); + } +} + +Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob( + const SimulationJob &fullModelSimulationJob) { + SimulationResults fullModelResults = + simulator.executeSimulation(fullModelSimulationJob, false, false); + fullModelResults.simulationLabel = "fullModel"; + if (draw) { + drawWorldAxes(); + fullModelResults.draw(fullModelSimulationJob); + } + computeDesiredReducedModelDisplacements(fullModelResults, + optimalReducedModelDisplacements); + computeReducedModelSimulationJob(fullModelSimulationJob, + reducedModelSimulationJob); + + // Set initial guess of solution + Eigen::VectorXd initialGuess(4); + auto propertiesOfFirstBeamOfFullModel = + pFullModelElementalMesh->elements[0].properties; + initialParameters(0) = propertiesOfFirstBeamOfFullModel.A; + initialParameters(1) = propertiesOfFirstBeamOfFullModel.J; + initialParameters(2) = propertiesOfFirstBeamOfFullModel.I2; + initialParameters(3) = propertiesOfFirstBeamOfFullModel.I3; + const double stifnessFactor = 2; + initialGuess(0) = stifnessFactor; + initialGuess(1) = stifnessFactor; + initialGuess(2) = stifnessFactor; + initialGuess(3) = stifnessFactor; + + const bool useGradientDescent = false; + if (useGradientDescent) { + + // gdc::GradientDescent, + // OptimizationCallback> + gdc::GradientDescent, + OptimizationCallback> + // gdc::GradientDescent, + // OptimizationCallback> + // gdc::GradientDescent, + // OptimizationCallback> + optimizer; + + // Turn verbosity on, so the optimizer prints status updates after each + // iteration. + optimizer.setVerbosity(1); + + // Set initial guess. + matplot::xlabel("Optimization iterations"); + matplot::ylabel("Objective value"); + // matplot::figure(false); + matplot::grid(matplot::on); + // Start the optimization + auto result = optimizer.minimize(initialGuess); + + std::cout << "Done! Converged: " << (result.converged ? "true" : "false") + << " Iterations: " << result.iterations << std::endl; + + // do something with final function value + std::cout << "Final fval: " << result.fval << std::endl; + + // do something with final x-value + std::cout << "Final xval: " << result.xval.transpose() << std::endl; + + SimulationResults reducedModelOptimizedResults = + simulator.executeSimulation(reducedModelSimulationJob); + reducedModelOptimizedResults.simulationLabel = "reducedModel"; + reducedModelOptimizedResults.draw(reducedModelSimulationJob); + + return result.xval; + } else { // use bobyqa + double (*pObjectiveFunction)(long, const double *) = &objective; + const size_t n = 4; + const size_t npt = 8; + assert(npt <= (n + 1) * (n + 2) / 2 && npt >= n + 2); + assert(npt < 2 * n + 1 && "The choice of the number of interpolation " + "conditions is not recommended."); + std::vector x{initialGuess(0), initialGuess(1), initialGuess(2), + initialGuess(3)}; + std::vector xLow(x.size(), -100); + std::vector xUpper(x.size(), 100); + const double maxX = *std::max_element( + x.begin(), x.end(), + [](const double &a, const double &b) { return abs(a) < abs(b); }); + const double rhobeg = std::min(0.95, 0.2 * maxX); + const double rhoend = rhobeg * 1e-6; + const size_t wSize = (npt + 5) * (npt + n) + 3 * n * (n + 5) / 2; + std::vector w(wSize); + bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(), + rhobeg, rhoend, 100, w.data()); + + std::cout << "Final objective value:" << objective(n, x.data()) + << std::endl; + + Eigen::VectorXd eigenX(x.size(), 1); + + for (size_t xi = 0; xi < x.size(); xi++) { + eigenX(xi) = x[xi]; + } + + SimulationResults reducedModelOptimizedResults = + simulator.executeSimulation(reducedModelSimulationJob); + reducedModelOptimizedResults.simulationLabel = "reducedModel"; + reducedModelOptimizedResults.draw(reducedModelSimulationJob); + + return eigenX; + } +} + +std::vector ReducedModelOptimizer::createScenarios( + const std::shared_ptr &pMesh) { + std::vector scenarios; + // Pull up + ////1Vertex + std::unordered_map> + pullUp1_FixedVertices; + pullUp1_FixedVertices[3] = std::unordered_set{0, 1, 2, 3, 4, 5}; + std::unordered_map pullUp1_NodalForces{ + {15, {0, 0, 100, 0, 0, 0}}}; + scenarios.push_back({pMesh, pullUp1_FixedVertices, pullUp1_NodalForces, {}}); + + // /// 2Vertices + // std::unordered_map> + // pullUp1_FixedVertices; + // saddleFixedVertices[3] = std::unordered_set{0, 1, 2}; + // saddleFixedVertices[7] = std::unordered_set{0, 1, 2}; + // std::unordered_map pullUp1_NodalForces{ + // {15, {0, 0, 0, 0, -900, 0}}, {3, {0, 0, 0, 0, 900, 0}}, + // {7, {0, 0, 0, 700, 0, 0}}, {11, {0, 0, 0, 700, 0, 0}}, + // {19, {0, 0, 0, -700, 0, 0}}, {23, {0, 0, 0, -700, 0, 0}}}; + // scenarios.push_back({pMesh, pullUp1_FixedVertices, pullUp1_NodalForces, + // {}}); + + /// 3Vertices + // std::unordered_map> + // pullUp1_FixedVertices; + // saddleFixedVertices[3] = std::unordered_set{0, 1, 2}; + // saddleFixedVertices[7] = std::unordered_set{0, 1, 2}; + // saddleFixedVertices[11] = std::unordered_set{0, 1, 2}; + // saddleFixedVertices[15] = std::unordered_set{0, 1, 2}; + // saddleFixedVertices[19] = std::unordered_set{0, 1, 2}; + // saddleFixedVertices[23] = std::unordered_set{0, 1, 2}; + // std::unordered_map pullUp1_NodalForces{ + // {15, {0, 0, 0, 0, -900, 0}}, {3, {0, 0, 0, 0, 900, 0}}, + // {7, {0, 0, 0, 700, 0, 0}}, {11, {0, 0, 0, 700, 0, 0}}, + // {19, {0, 0, 0, -700, 0, 0}}, {23, {0, 0, 0, -700, 0, 0}}}; + // scenarios.push_back({pMesh, pullUp1_FixedVertices, pullUp1_NodalForces, + // {}}); + + // Single moment + // std::unordered_map> + // fixedVertices; fixedVertices[3] = std::unordered_set{0, 1, 2, 3, + // 4, 5}; fixedVertices[15] = std::unordered_set{0, 1, 2}; + // std::unordered_map nodalForces{ + // {15, {0, 0, 0, 0, 700, 0}}}; + // scenarios.push_back({pMesh, saddleFixedVertices, saddleNodalForces, {}}); + // scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); + + // Saddle + // std::unordered_map> + // saddle_fixedVertices; + // // saddle_fixedVertices[3] = std::unordered_set{0, 1, 2}; + // saddle_fixedVertices[7] = std::unordered_set{0, 1, 2}; + // saddle_fixedVertices[11] = std::unordered_set{0, 1, 2}; + // // saddle_fixedVertices[15] = std::unordered_set{0, 1, 2}; + // // saddle_fixedVertices[19] = std::unordered_set{0, 1, 2}; + // // saddle_fixedVertices[23] = std::unordered_set{0, 1, 2}; + // std::unordered_map saddle_nodalForces{ + // {15, {0, 0, 0, 0, -1400, 0}}, {3, {0, 0, 0, 0, 1400, 0}}, + // {7, {0, 0, 0, 700, 0, 0}}, {11, {0, 0, 0, 700, 0, 0}}, + // {19, {0, 0, 0, -700, 0, 0}}, {23, {0, 0, 0, -700, 0, 0}}}; + // // scenarios.push_back({pMesh, saddleFixedVertices, saddleNodalForces, + // {}}); std::unordered_map + // saddle_forcedDisplacements; scenarios.push_back({pMesh, + // saddle_fixedVertices, saddle_nodalForces, + // saddle_forcedDisplacements}); + + return scenarios; +} + +Eigen::VectorXd ReducedModelOptimizer::optimize() { + std::vector simulationJobs = + createScenarios(pFullModelElementalMesh); + std::vector results; + for (const SimulationJob &job : simulationJobs) { + // fullModelSimulationJob.draw(); + auto result = optimizeForSimulationJob(job); + results.push_back(result); + } + + if (results.empty()) { + return Eigen::VectorXd(); + } + + return results[0]; +} diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp new file mode 100644 index 0000000..675d87f --- /dev/null +++ b/src/reducedmodeloptimizer.hpp @@ -0,0 +1,57 @@ +#ifndef REDUCEDMODELOPTIMIZER_HPP +#define REDUCEDMODELOPTIMIZER_HPP + +#include "beamformfinder.hpp" +#include "edgemesh.hpp" +#include "elementalmesh.hpp" +#include + +using FullModelVertexIndex = VertexIndex; +using ReducedModelVertexIndex = VertexIndex; + +struct SimulationScenarios { + std::vector scenarios; + std::shared_ptr pMesh; + SimulationScenarios(std::shared_ptr pMesh) : pMesh(pMesh) {} + + std::vector getSimulationJobs() const { return scenarios; } +}; + +class ReducedModelOptimizer { + std::shared_ptr reducedModelElementalMesh; + std::shared_ptr pFullModelElementalMesh; + const std::unordered_map + &reducedToFullViMap; + std::unordered_map + fullToReducedViMap; + std::unordered_map + oppositeVertexMap; + +public: + Eigen::VectorXd optimize(); + double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const; + + ReducedModelOptimizer( + ConstVCGEdgeMesh &fullModel, ConstVCGEdgeMesh &reducedModel, + const std::unordered_map + &viMap, + const std::vector> &oppositeVertices); + void computeReducedModelSimulationJob( + const SimulationJob &simulationJobOfFullModel, + SimulationJob &simulationJobOfReducedModel); + + SimulationJob + getReducedSimulationJob(const SimulationJob &fullModelSimulationJob); + +private: + void computeDesiredReducedModelDisplacements( + const SimulationResults &fullModelResults, + Eigen::MatrixX3d &optimaldisplacementsOfReducedModel); + Eigen::VectorXd + optimizeForSimulationJob(const SimulationJob &fullModelSimulationJob); + std::vector + createScenarios(const std::shared_ptr &pMesh); +}; + +#endif // REDUCEDMODELOPTIMIZER_HPP