diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..e69de29 diff --git a/.travis.yml b/.travis.yml index c0c1fdc..6192bd7 100644 --- a/.travis.yml +++ b/.travis.yml @@ -53,7 +53,12 @@ install: script: - source activate testenv - nosetests --with-coverage --cover-package=prox_tv -- bash -x ./.travis/deploy.sh +# - bash -x ./.travis/deploy.sh +# Build C lib with CMake +- mkdir -p build && cd build +- cmake -DENABLE_TESTING:BOOL=ON ../ +- cmake --build . +- ctest -VV after_success: - source activate testenv - coveralls diff --git a/.travis/installconda.sh b/.travis/installconda.sh index dcb9ede..20a9314 100644 --- a/.travis/installconda.sh +++ b/.travis/installconda.sh @@ -22,6 +22,7 @@ export PATH="$HOME/miniconda3/bin:$PATH" rm miniconda.sh conda config --set always_yes yes --set changeps1 no conda update -q conda +pip install cmake # Create environment with specific python version conda create -n testenv python=${TRAVIS_PYTHON_VERSION} diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..f1f396c --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,147 @@ +cmake_minimum_required(VERSION 3.12) + +project(proxTV + VERSION 3.2.1 + DESCRIPTION "Toolbox for fast Total Variation proximity operators" + ) +message(STATUS "Configuring ${PROJECT_NAME}") +message(STATUS " version: ${proxTV_VERSION}") +message(STATUS " description: ${proxTV_DESCRIPTION}") + +# Update CMake module path to lookup proxTV custom CMake modules +list(INSERT CMAKE_MODULE_PATH 0 ${CMAKE_CURRENT_SOURCE_DIR}/cmake) + +# Add -fPIC +set(CMAKE_POSITION_INDEPENDENT_CODE ON) + +#------------------------------------------------------------------------------ +if(NOT DEFINED proxTV_INSTALL_DEVELOPMENT) + option(proxTV_INSTALL_DEVELOPMENT "Install development files" ON) +endif() +if(NOT DEFINED proxTV_ENABLE_TESTING) + option(proxTV_ENABLE_TESTING "Compile tests" OFF) +endif() +if(NOT DEFINED proxTV_USE_LAPACK) + option(proxTV_USE_LAPACK "Use LAPACK and LAPACKE instead of EIGEN" ON) +endif() + +#------------------------------------------------------------------------------ +# Set a default build type if none was specified +if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) + message(STATUS "Setting build type to 'Release' as none was specified.") + set(CMAKE_BUILD_TYPE "Release" CACHE STRING "Choose the type of build." FORCE) + mark_as_advanced(CMAKE_BUILD_TYPE) + # Set the possible values of build type for cmake-gui + set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release" "MinSizeRel" "RelWithDebInfo") +endif() + +#------------------------------------------------------------------------------ +# Install directories +if(NOT DEFINED proxTV_INSTALL_BIN_DIR) + set(proxTV_INSTALL_BIN_DIR bin) +endif() +if(NOT DEFINED proxTV_INSTALL_LIB_DIR) + set(proxTV_INSTALL_LIB_DIR lib) +endif() +if(NOT DEFINED proxTV_INSTALL_INCLUDE_DIR) + set(proxTV_INSTALL_INCLUDE_DIR include) +endif() +if(NOT DEFINED proxTV_INSTALL_CMAKE_DIR) + set(proxTV_INSTALL_CMAKE_DIR ${proxTV_INSTALL_LIB_DIR}/cmake/proxTV) +endif() + +message(STATUS "proxTV_INSTALL_INCLUDE_DIR: ${proxTV_INSTALL_INCLUDE_DIR}") +message(STATUS "proxTV_INSTALL_LIB_DIR: ${proxTV_INSTALL_LIB_DIR}") +message(STATUS "proxTV_INSTALL_CMAKE_DIR: ${proxTV_INSTALL_CMAKE_DIR}") +#------------------------------------------------------------------------------ +# External dependencies + +if(proxTV_USE_LAPACK) + find_package(LAPACKE REQUIRED) + message(STATUS "Lapacke libraries: ${LAPACKE_LIBRARIES}") + + find_package(LAPACK REQUIRED) + message(STATUS "Lapack libraries: ${LAPACK_LIBRARIES}") +else() + # No need to find Eigen3 if the Eigen3 target exists already. + # This happens when a project is fetching Eigen before fetching proxTV with add_subdirectory + # Particularly, this allows to build the module inside ITK when Module_TotalVariation is ON + if(NOT TARGET Eigen3::Eigen3) + find_package(Eigen3 REQUIRED) + endif() + get_target_property(EIGEN_INCLUDE_DIR Eigen3::Eigen + INTERFACE_INCLUDE_DIRECTORIES) + message(STATUS "Eigen Found: ${EIGEN_INCLUDE_DIR}") +endif() + +set(THREADS_PREFER_PTHREAD_FLAG 1) +find_package(Threads) + +find_package(OpenMP) +message(STATUS "OpenMP found: ${OpenMP_FOUND}") + +#------------------------------------------------------------------------------ +# Add libraries + +add_subdirectory(src) + +#------------------------------------------------------------------------------ +# Testing +if(proxTV_ENABLE_TESTING) + enable_testing() + add_subdirectory(test) +endif() + +#------------------------------------------------------------------------------ +# Configure proxTVConfigVersion.cmake common to build and install tree +include(CMakePackageConfigHelpers) +set(config_version_file "${proxTV_BINARY_DIR}/proxTVConfigVersion.cmake") +write_basic_package_version_file(${config_version_file} + VERSION ${proxTV_VERSION} + COMPATIBILITY SameMajorVersion + ) + +#------------------------------------------------------------------------------ +# Export 'proxTVTargets.cmake' for a build tree +export(TARGETS proxTV + FILE ${PROJECT_BINARY_DIR}/proxTVTargets.cmake + NAMESPACE proxTV:: + ) + +# Configure 'proxTVConfig.cmake' for a build tree +include(CMakePackageConfigHelpers) +set(build_config ${PROJECT_BINARY_DIR}/proxTVConfig.cmake) +configure_package_config_file( + cmake/proxTVConfig.cmake.in + ${build_config} + INSTALL_DESTINATION ${PROJECT_BINARY_DIR} + INSTALL_PREFIX ${PROJECT_BINARY_DIR} + NO_CHECK_REQUIRED_COMPONENTS_MACRO + ) + +#------------------------------------------------------------------------------ +# Configure 'proxTVConfig.cmake' for an install tree +if(proxTV_INSTALL_DEVELOPMENT) + set(install_config ${PROJECT_BINARY_DIR}/CMakeFiles/proxTVConfig.cmake) + configure_package_config_file( + cmake/proxTVConfig.cmake.in + ${install_config} + INSTALL_DESTINATION ${proxTV_INSTALL_CMAKE_DIR} + NO_CHECK_REQUIRED_COMPONENTS_MACRO + ) + + # Install 'proxTVTargets.cmake' + install(EXPORT proxTVTargets + FILE proxTVTargets.cmake + NAMESPACE proxTV:: + DESTINATION ${proxTV_INSTALL_CMAKE_DIR} + COMPONENT Development + ) + + # Install config files + install( + FILES ${config_version_file} ${install_config} + DESTINATION ${proxTV_INSTALL_CMAKE_DIR} + COMPONENT Development + ) +endif() diff --git a/README.md b/README.md index 6ba6280..dbeca13 100644 --- a/README.md +++ b/README.md @@ -36,6 +36,33 @@ More technically, the library provides efficient solvers for the following Total | Anisotropic Total Variation on a 3-dimensional signal (video denoising) | ![alt tag](docs/img/TV3D.png) | | Generalized N-dimensional Anisotropic Total Variation (tensor denoising) | ![alt tag](docs/img/TVND.png), with X(di) every possible 1-dimensional slice of X following dimension di.| +## C interface + +You can generate a **c** static or dynamic library using **cmake**. If `libproxTV` is not provided by your package-manager, install it from source: + + mkdir proxTV-dev ; cd proxTV-dev + git clone https://github.com/albarji/proxTV proxTV + mkdir build ; cd build + cmake ../proxTV + make -j4 + make install + +The required dependencies are `lapack` and `lapacke`, the c-interface for `lapack`, and optionally, but recommended: `OpenMP` with `pthreads` for multi-threading support. + +You can provide extra options to `cmake` via the command line or a gui (i.e `ccmake`). + + cmake ../proxTV -DBUILD_SHARED_LIBS:BOOL=ON -DCMAKE_INSTALL_PREFIX=/opt/ -DENABLE_TESTING:BOOL=ON + +To use proxTV in your `cmake` project just write in your `CMakeLists.txt`: + + find_package(proxTV) + add_executable(foo main.cpp) + target_link_libraries(foo PUBLIC proxTV::proxTV) + +That will propagate all the dependencies of proxTV to your target. If you haven't installed proxTV in a system folder, you have to point to the installation directory when configuring your project with `cmake`. + + cmake /path/my_project_source_folder -DproxTV_DIR:PATH="/proxTV_install_folder/lib/cmake/proxTV" + ## Python interface ### Install diff --git a/cmake/FindLAPACKE.cmake b/cmake/FindLAPACKE.cmake new file mode 100644 index 0000000..b0e9639 --- /dev/null +++ b/cmake/FindLAPACKE.cmake @@ -0,0 +1,190 @@ +#.rst: +# FindLAPACKE +# ------------- +# +# Find the LAPACKE library +# +# Using LAPACKE: +# +# :: +# +# find_package(LAPACKE REQUIRED) +# include_directories(${LAPACKE_INCLUDE_DIRS}) +# add_executable(foo foo.cc) +# target_link_libraries(foo ${LAPACKE_LIBRARIES}) +# +# This module sets the following variables: +# +# :: +# +# LAPACKE_FOUND - set to true if the library is found +# LAPACKE_INCLUDE_DIRS - list of required include directories +# LAPACKE_LIBRARIES - list of libraries to be linked +# LAPACKE_VERSION_MAJOR - major version number +# LAPACKE_VERSION_MINOR - minor version number +# LAPACKE_VERSION_PATCH - patch version number +# LAPACKE_VERSION_STRING - version number as a string (ex: "0.2.18") + +#============================================================================= +# Copyright 2016 Hans J. Johnson +# +# Distributed under the OSI-approved BSD License (the "License") +# +# This software is distributed WITHOUT ANY WARRANTY; without even the +# implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +#============================================================================= +# +set(LAPACKE_SEARCH_PATHS + ${LAPACKE_DIR} + $ENV{LAPACKE_DIR} + $ENV{CMAKE_PREFIX_PATH} + ${CMAKE_PREFIX_PATH} + /usr + /usr/local + /usr/local/opt/lapack ## Mac Homebrew install path + /opt/LAPACKE +) +message(STATUS "LAPACKE_SEARCH_PATHS: ${LAPACKE_SEARCH_PATHS}") + +set(CMAKE_PREFIX_PATH ${LAPACKE_SEARCH_PATHS}) +list(REMOVE_DUPLICATES CMAKE_PREFIX_PATH) + +## First try to find LAPACKE with NO_MODULE, +## As of 20160706 version 0.2.18 there is limited cmake support for LAPACKE +## that is not as complete as this version, if found, use it +## to identify the LAPACKE_VERSION_STRING and improve searching. +find_package(LAPACKE NO_MODULE QUIET) +if(LAPACKE_FOUND) + if(EXISTS ${LAPACKE_DIR}/lapacke-config-version.cmake) + include(${LAPACKE_DIR}/lapacke-config-version.cmake) + set(LAPACKE_VERSION_STRING ${PACKAGE_VERSION}) + unset(PACKAGE_VERSION) # Use cmake conventional naming + endif() + find_package(LAPACK NO_MODULE QUIET) #Require matching versions here! + find_package(BLAS NO_MODULE QUIET) #Require matching versions here! +endif() + +################################################################################################## +### First search for headers +find_path(LAPACKE_CBLAS_INCLUDE_DIR + NAMES cblas.h + PATHS ${LAPACKE_SEARCH_PATHS} + PATH_SUFFIXES include include/lapack) +find_path(LAPACKE_LAPACKE_INCLUDE_DIR + NAMES lapacke.h + PATHS ${LAPACKE_SEARCH_PATHS} + PATH_SUFFIXES include) + +################################################################################################## +### Second, search for libraries +set(PATH_SUFFIXES_LIST + lib64 + lib +) +find_library(LAPACKE_LIB + NAMES lapacke + PATHS ${LAPACKE_SEARCH_PATHS} + PATH_SUFFIXES ${PATH_SUFFIXES_LIST}) +find_library(CBLAS_LIB + NAMES cblas + PATHS ${LAPACKE_SEARCH_PATHS} + PATH_SUFFIXES ${PATH_SUFFIXES_LIST}) +find_library(LAPACK_LIB + NAMES lapack + PATHS ${LAPACKE_SEARCH_PATHS} + PATH_SUFFIXES ${PATH_SUFFIXES_LIST}) +find_library(BLAS_LIB + NAMES blas + PATHS ${LAPACKE_SEARCH_PATHS} + PATH_SUFFIXES ${PATH_SUFFIXES_LIST}) + +## TODO: Get version components +# ------------------------------------------------------------------------ +# Extract version information +# ------------------------------------------------------------------------ + +# WARNING: We may not be able to determine the version of some LAPACKE +set(LAPACKE_VERSION_MAJOR 0) +set(LAPACKE_VERSION_MINOR 0) +set(LAPACKE_VERSION_PATCH 0) +if(LAPACKE_VERSION_STRING) + string(REGEX REPLACE "([0-9]+).([0-9]+).([0-9]+)" "\\1" LAPACKE_VERSION_MAJOR "${LAPACKE_VERSION_STRING}") + string(REGEX REPLACE "([0-9]+).([0-9]+).([0-9]+)" "\\2" LAPACKE_VERSION_MINOR "${LAPACKE_VERSION_STRING}") + string(REGEX REPLACE "([0-9]+).([0-9]+).([0-9]+)" "\\3" LAPACKE_VERSION_PATCH "${LAPACKE_VERSION_STRING}") +endif() + +#====================== +# Checks 'REQUIRED', 'QUIET' and versions. +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(LAPACKE FOUND_VAR LAPACKE_FOUND + REQUIRED_VARS LAPACKE_CBLAS_INCLUDE_DIR + LAPACKE_LAPACKE_INCLUDE_DIR + LAPACKE_LIB + LAPACK_LIB + CBLAS_LIB + BLAS_LIB + VERSION_VAR LAPACKE_VERSION_STRING +) + +if (LAPACKE_FOUND) + set(LAPACKE_INCLUDE_DIRS ${LAPACKE_CBLAS_INCLUDE_DIR} ${LAPACKE_CBLAS_INCLUDE_DIR}) + list(REMOVE_DUPLICATES LAPACKE_INCLUDE_DIRS) + if("${CMAKE_C_COMPILER_ID}" MATCHES ".*Clang.*" OR + "${CMAKE_C_COMPILER_ID}" MATCHES ".*GNU.*" OR + "${CMAKE_C_COMPILER_ID}" MATCHES ".*Intel.*" + ) #NOT MSVC + set(MATH_LIB m) + endif() + list(APPEND LAPACKE_LIBRARIES ${LAPACKE_LIB} ${LAPACK_LIB} ${BLAS_LIB} ${CBLAS_LIB}) + # Check for a common combination, and find required gfortran support libraries + + if(1) + if("${CMAKE_C_COMPILER_ID}" MATCHES ".*Clang.*" AND "${CMAKE_Fortran_COMPILER_ID}" MATCHES "GNU") + message(STATUS "\n\n WARNING: ${CMAKE_C_COMPILER} identified as ${CMAKE_C_COMPILER_ID}\n" + "AND: ${CMAKE_Fortran_COMPILER} identified as ${CMAKE_Fortran_COMPILER_ID}\n" + "\n" + "may be require special configurations. The most common is the need to" + "explicitly link C programs against the gfortran support library.") + + endif() + else() + ## This code automated code is hard to determine if it is robust in many different environments. + # Check for a common combination, and find required gfortran support libraries + if("${CMAKE_C_COMPILER_ID}" MATCHES ".*Clang.*" AND "${CMAKE_Fortran_COMPILER_ID}" MATCHES "GNU") + include(FortranCInterface) + FortranCInterface_VERIFY() + if(NOT FortranCInterface_VERIFIED_C) + message(FATAL_ERROR "C and fortran compilers are not compatible:\n${CMAKE_Fortran_COMPILER}:${CMAKE_C_COMPILER}") + endif() + + execute_process(COMMAND ${CMAKE_Fortran_COMPILER} -print-file-name=libgfortran.a OUTPUT_VARIABLE FORTRANSUPPORTLIB ERROR_QUIET) + string(STRIP ${FORTRANSUPPORTLIB} FORTRANSUPPORTLIB) + if(EXISTS "${FORTRANSUPPORTLIB}") + list(APPEND LAPACKE_LIBRARIES ${FORTRANSUPPORTLIB}) + message(STATUS "Appending fortran support lib: ${FORTRANSUPPORTLIB}") + else() + message(FATAL_ERROR "COULD NOT FIND libgfortran.a support library:${FORTRANSUPPORTLIB}:") + endif() + endif() + endif() + list(APPEND LAPACKE_LIBRARIES ${MATH_LIB}) +endif() + +mark_as_advanced( + LAPACKE_FOUND + LAPACKE_INCLUDE_DIRS + LAPACKE_LIBRARIES + LAPACKE_VERSION_MAJOR + LAPACKE_VERSION_MINOR + LAPACKE_VERSION_PATCH + LAPACKE_VERSION_STRING +) + +## For debugging +message(STATUS "LAPACKE_FOUND :${LAPACKE_FOUND}: - set to true if the library is found") +message(STATUS "LAPACKE_INCLUDE_DIRS :${LAPACKE_INCLUDE_DIRS}: - list of required include directories") +message(STATUS "LAPACKE_LIBRARIES :${LAPACKE_LIBRARIES}: - list of libraries to be linked") +message(STATUS "LAPACKE_VERSION_MAJOR :${LAPACKE_VERSION_MAJOR}: - major version number") +message(STATUS "LAPACKE_VERSION_MINOR :${LAPACKE_VERSION_MINOR}: - minor version number") +message(STATUS "LAPACKE_VERSION_PATCH :${LAPACKE_VERSION_PATCH}: - patch version number") +message(STATUS "LAPACKE_VERSION_STRING :${LAPACKE_VERSION_STRING}: - version number as a string") diff --git a/cmake/README.md b/cmake/README.md new file mode 100644 index 0000000..be07c7c --- /dev/null +++ b/cmake/README.md @@ -0,0 +1,7 @@ +Modern lapacke provides a lapacke-config.cmake (minimum version unconfirmed, but 3.7.1 does) +so FindLAPACKE.cmake wouldn't be neccesary, (see #38) +however this isn't reliable for older/other lapacke versions, so we provide a FindModule from: https://github.com/mrirecon/bart/blob/master/cmake/FindLAPACKE.cmake. + +Download with `curl -O https://raw.githubusercontent.com/mrirecon/bart/master/cmake/FindLAPACKE.cmake` + +This FindLapacke handles properly the case a lapacke-config.cmake exists in the system. diff --git a/cmake/proxTVConfig.cmake.in b/cmake/proxTVConfig.cmake.in new file mode 100644 index 0000000..2a641a7 --- /dev/null +++ b/cmake/proxTVConfig.cmake.in @@ -0,0 +1,17 @@ +@PACKAGE_INIT@ + +set_and_check(proxTV_TARGETS "${CMAKE_CURRENT_LIST_DIR}/proxTVTargets.cmake") + +if(NOT TARGET proxTV::proxTV) + include(${proxTV_TARGETS}) +endif() + +include(CMakeFindDependencyMacro) +if(@proxTV_USE_LAPACK@) + find_dependency(LAPACKE) + find_dependency(LAPACK) +else() + find_dependency(Eigen3) +endif() +find_dependency(OpenMP) +find_dependency(Threads) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt new file mode 100644 index 0000000..f320421 --- /dev/null +++ b/src/CMakeLists.txt @@ -0,0 +1,99 @@ +set(headers + condat_fast_tv.h + general.h + johnsonRyanTV.h + lapackFunctionsWrap.h + LPopt.h + TVmacros.h + TVopt.h + utils.h + ) + +set(sources + condat_fast_tv.cpp + johnsonRyanTV.cpp + lapackFunctionsWrap.cpp + LPopt.cpp + TV2Dopt.cpp + TV2DWopt.cpp + TVgenopt.cpp + TVL1opt.cpp + TVL1opt_hybridtautstring.cpp + TVL1opt_kolmogorov.cpp + TVL1opt_tautstring.cpp + TVL1Wopt.cpp + TVL2opt.cpp + TVLPopt.cpp + TVNDopt.cpp + utils.cpp + ) + +add_library(proxTV ${sources}) +set_target_properties(proxTV PROPERTIES WINDOWS_EXPORT_ALL_SYMBOLS 1) +target_compile_definitions(proxTV PUBLIC NOMATLAB) +if(proxTV_USE_LAPACK) + target_compile_definitions(proxTV PUBLIC PROXTV_USE_LAPACK) +endif() + +# Create directory with headers at build and install time. +file(COPY ${headers} + DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/include") +install(FILES ${headers} + DESTINATION ${proxTV_INSTALL_INCLUDE_DIR}) + +target_include_directories(proxTV PUBLIC + $ + $ + $ + ) + +# Pthreads +target_link_libraries(proxTV PUBLIC + ${CMAKE_THREAD_LIBS_INIT} + ) + +# Lapacke, Lapack +if(proxTV_USE_LAPACK) + target_include_directories(proxTV PUBLIC + ${LAPACKE_INCLUDE_DIRS} + ) + target_link_libraries(proxTV PUBLIC + ${LAPACK_LIBRARIES} + ${LAPACKE_LIBRARIES} + ) +else() + target_link_libraries(proxTV PUBLIC Eigen3::Eigen) +endif() + +# OpenMP +if(OpenMP_CXX_FOUND AND NOT MSVC) + # target_compile_definitions(proxTV PRIVATE _OPENMP) + target_link_libraries(proxTV PRIVATE OpenMP::OpenMP_CXX) +endif() + +install(TARGETS proxTV EXPORT proxTVTargets + LIBRARY DESTINATION ${proxTV_INSTALL_LIB_DIR} + ARCHIVE DESTINATION ${proxTV_INSTALL_LIB_DIR} + RUNTIME DESTINATION ${proxTV_INSTALL_BIN_DIR} + INCLUDES DESTINATION ${proxTV_INSTALL_INCLUDE_DIR} + ) + +# This saves as to develop a hacky FindproxTV for other to use the library. +# A regular find_package(proxTV) will look for these targets. +# It needs an extra proxTVConfig.cmake to handle dependencies (provided in cmake directory) +install(EXPORT proxTVTargets + FILE proxTVTargets.cmake + NAMESPACE proxTV:: + DESTINATION ${proxTV_INSTALL_CMAKE_DIR} + ) + +# Generate and install proxTVConfigVersion.cmake +include(CMakePackageConfigHelpers) +write_basic_package_version_file("proxTVConfigVersion.cmake" + VERSION ${proxTV_VERSION} + COMPATIBILITY SameMajorVersion + ) +install(FILES + ${CMAKE_CURRENT_BINARY_DIR}/proxTVConfigVersion.cmake + DESTINATION ${proxTV_INSTALL_CMAKE_DIR} + ) diff --git a/src/LPopt.cpp b/src/LPopt.cpp index 4203f1b..e70d981 100644 --- a/src/LPopt.cpp +++ b/src/LPopt.cpp @@ -211,7 +211,10 @@ int PN_LPinf(double *y,double lambda,double *x,double *info,int n,Workspace *ws) */ int PN_LPp(double *y,double lambda,double *x,double *info,int n,double p,Workspace *ws,int positive,double objGap){ double *g=NULL,*d=NULL,*xnorm=NULL,*auxv=NULL; - double stop,stop2,q,nx,f,fupdate,aux,c,den,xp1vGrad,gRd,delta,prevDelta,improve,rhs,grad0,gap,epsilon; + double stop,stop2,q,nx,f,aux,c,den,xp1vGrad,gRd,delta,prevDelta,rhs,gap,epsilon; + double fupdate = 0; + double improve = 0; + double grad0 = 0; int *inactive=NULL,*signs=NULL; int i,j,iters,recomp,found,nI; short updateKind; @@ -745,7 +748,7 @@ int PN_LPp(double *y,double lambda,double *x,double *info,int n,double p,Workspa @argument lambda multiplier of the Lp norm @argument norm precomputed Lp norm of x **/ -double PN_LPpGap(double *x, double *y, double *diff, int n, double q, double lambda, double norm) { +double PN_LPpGap(double *, double *y, double *diff, int n, double q, double lambda, double norm) { /* Compute dual norm */ double dualnorm = LPnorm(diff,n,q); @@ -1016,7 +1019,7 @@ void solveLinearLP(double *z, int n, double p, double lambda, double *s) { #endif // The solution approximately lies at the opposite corner of the l1 ball double largest = 0, val; - int ilargest; + int ilargest = 0; // Find largest entry of z, in absolute value for ( i = 0 ; i < n ; i++ ) { diff --git a/src/TV2DWopt.cpp b/src/TV2DWopt.cpp index 76cb5af..b860cfc 100644 --- a/src/TV2DWopt.cpp +++ b/src/TV2DWopt.cpp @@ -46,7 +46,7 @@ void DR_rowsPass(size_t M, size_t N, double* input, double* output, double* ref, int DR2L1W_TV(size_t M, size_t N, double*unary, double*W1, double*W2, double*s, int nThreads, int maxit, double* info) { - int i; + size_t i; double *t = NULL; double *tb = NULL; Workspace **ws = NULL; @@ -67,7 +67,7 @@ int DR2L1W_TV(size_t M, size_t N, double*unary, double*W1, double*W2, double*s, if (nThreads < 1) nThreads = 1; omp_set_num_threads(nThreads); - maxDim = (M > N) ? M : N; + maxDim = int( (M > N) ? M : N ); // Alloc memory for algorithm */ t = (double*) malloc(sizeof(double)*M*N); @@ -152,7 +152,6 @@ int DR2L1W_TV(size_t M, size_t N, double*unary, double*W1, double*W2, double*s, void DR_columnsPass(size_t M, size_t N, double* input, double* output, double* W, Workspace **ws) { #pragma omp parallel shared(M,N,input,output,W,ws) default(none) { - int i,j; // Get thread number int id = omp_get_thread_num(); // Get corresponding workspace @@ -161,7 +160,7 @@ void DR_columnsPass(size_t M, size_t N, double* input, double* output, double* W // Run 1-d solvers in parallel on each column of the input #pragma omp for - for (j=0; j < N; j++) { + for (int j=0; j < N; j++) { resetWorkspace(wsi); // Array for weights double* wline = getDoubleWorkspace(wsi); @@ -205,7 +204,7 @@ void DR_rowsPass(size_t M, size_t N, double* input, double* output, double* ref, // Array for weights double* wline = getDoubleWorkspace(wsi); // Prepare weights - int idx; + size_t idx; for ( idx = j, i = 0 ; i < N-1 ; i++, idx+=M ) wline[i] = W[idx]; // Prepare inputs, considering displacement from reference signal @@ -229,13 +228,12 @@ void DR_rowsPass(size_t M, size_t N, double* input, double* output, double* ref, @param W weights of the TV regularization @param ws Workspace to use for the computation */ -void DR_proxDiff(size_t n, double* input, double* output, double* W, Workspace *ws) { - int i; +void DR_proxDiff(size_t n, double* input, double* output, double* W, Workspace *) { // Compute proximity - tautString_TV1_Weighted(input, W, output, n); + tautString_TV1_Weighted(input, W, output, (int) n); // Return differences between input and proximity output - for (i=0; i < n; i++) + for (size_t i=0; i < n; i++) output[i] = input[i] - output[i]; } diff --git a/src/TV2Dopt.cpp b/src/TV2Dopt.cpp index a89045f..0c5ccfb 100644 --- a/src/TV2Dopt.cpp +++ b/src/TV2Dopt.cpp @@ -168,9 +168,9 @@ int PD2_TV(double *y,double *lambdas,double *norms,double *dims,double *x,double #ifdef DEBUG fprintf(DEBUG_FILE,"··········Penalty 0··········\n"); #endif - d = dims[0]-1; + d = int(dims[0]-1); /* Run 1-dimensional prox operator over each 1-dimensional slice along the specified dimension (parallelized) */ - #pragma omp parallel shared(ws,nSlices,ns,d,incs,x,p,lambdas,z,norms,DEBUG_FILE) private(j,k,idx1,idx2) default(none) + #pragma omp parallel shared(ws,nSlices,ns,d,incs,x,p,lambdas,z,norms) private(j,k,idx1,idx2) default(none) { /* Get thread number */ int id = omp_get_thread_num(); @@ -187,16 +187,6 @@ int PD2_TV(double *y,double *lambdas,double *norms,double *dims,double *x,double for(k=0,idx2=0 ; kin[k] = x[idx1+idx2]+p[idx1+idx2]; - #ifdef DEBUG - { - int dbgi; - fprintf(DEBUG_FILE,"Slice %d: ",j); - for(dbgi=0;dbgiin[dbgi]); - fprintf(DEBUG_FILE,"\n"); - } - #endif - /* Apply 1-dimensional solver */ resetWorkspace(wsi); TV(wsi->in, lambdas[0], wsi->out, NULL, ns[d], norms[0], wsi); @@ -214,13 +204,10 @@ int PD2_TV(double *y,double *lambdas,double *norms,double *dims,double *x,double /* Prox step for the second penalty term (if any) */ if(npen >= 2){ - #ifdef DEBUG - fprintf(DEBUG_FILE,"··········Penalty 1··········\n"); - #endif - d = dims[1]-1; + d = int(dims[1]-1); /* Run 1-dimensional prox operator over each 1-dimensional slice along the specified dimension (parallelized) */ - #pragma omp parallel shared(ws,nSlices,ns,d,incs,x,q,lambdas,z,norms,DEBUG_FILE) private(j,k,idx1,idx2) default(none) + #pragma omp parallel shared(ws,nSlices,ns,d,incs,x,q,lambdas,z,norms) private(j,k,idx1,idx2) default(none) { /* Get thread number */ int id = omp_get_thread_num(); @@ -237,16 +224,6 @@ int PD2_TV(double *y,double *lambdas,double *norms,double *dims,double *x,double for(k=0,idx2=0 ; kin[k] = z[idx1+idx2] + q[idx1+idx2]; - #ifdef DEBUG - { - int dbgi; - fprintf(DEBUG_FILE,"Slice %d: ",j); - for(dbgi=0;dbgiin[dbgi]); - fprintf(DEBUG_FILE,"\n"); - } - #endif - /* Apply 1-dimensional solver */ resetWorkspace(wsi); TV(wsi->in, lambdas[1], wsi->out, NULL, ns[d], norms[1], wsi); @@ -287,9 +264,6 @@ int PD2_TV(double *y,double *lambdas,double *norms,double *dims,double *x,double /* Termination check */ if(iters >= MAX_ITERS_PD){ - #ifdef DEBUG - fprintf(DEBUG_FILE,"(PD2_TV) WARNING: maximum number of iterations reached (%d).\n",MAX_ITERS_PD); - #endif if(info) info[INFO_RC] = RC_ITERS; } else if(info) info[INFO_RC] = RC_OK; @@ -353,12 +327,10 @@ int DR2_TV(size_t M, size_t N, double*unary, double W1, double W2, double norm1, double norm2, double*s, int nThreads, int maxit, double* info) { - int i; - double *ytr = NULL; + size_t i; double *t = NULL; double *tb = NULL; Workspace **ws = NULL; - int maxDim; #define FREE \ if(t) free(t); \ @@ -373,7 +345,7 @@ int DR2_TV(size_t M, size_t N, double*unary, double W1, double W2, if (nThreads < 1) nThreads = 1; omp_set_num_threads(nThreads); - maxDim = (M > N) ? M : N; + int maxDim = int( (M > N) ? M : N ); // Alloc memory for algorithm */ t = (double*) malloc(sizeof(double)*M*N); @@ -459,7 +431,6 @@ int DR2_TV(size_t M, size_t N, double*unary, double W1, double W2, void DR_columnsPass(size_t M, size_t N, double* input, double* output, double W, double norm, Workspace **ws) { #pragma omp parallel shared(M,N,input,output,W,norm,ws) default(none) { - int j; // Get thread number int id = omp_get_thread_num(); // Get corresponding workspace @@ -468,7 +439,7 @@ void DR_columnsPass(size_t M, size_t N, double* input, double* output, double W, // Run 1-d solvers in parallel on each column of the input #pragma omp for - for (j=0; j < N; j++) { + for (int j=0; j < N; j++) { resetWorkspace(wsi); // Prepare inputs memcpy(wsi->in, input+(M*j), sizeof(double)*M); @@ -510,7 +481,7 @@ void DR_rowsPass(size_t M, size_t N, double* input, double* output, double* ref, for (j=0; j < M; j++) { resetWorkspace(wsi); // Prepare inputs, considering displacement from reference signal - int idx; + size_t idx; for ( idx = j, i = 0 ; i < N ; i++, idx+=M ) wsi->in[i] = ref[idx]-input[idx]; // Compute prox difference for this row @@ -536,13 +507,12 @@ difference between input and output of prox. @param norm degree of TV @param ws Workspace to use for the computation */ -void DR_proxDiff(size_t n, double* input, double* output, double W, double norm, Workspace *ws) { - int i; +void DR_proxDiff(size_t n, double* input, double* output, double W, double norm, Workspace *) { // Compute proximity - TV(input, W, output, NULL, n , norm, NULL); + TV(input, W, output, NULL, (int) n , norm, NULL); // Return differences between input and proximity output - for (i=0; i < n; i++) + for (size_t i=0; i < n; i++) output[i] = input[i] - output[i]; } @@ -585,9 +555,10 @@ void DR_proxDiff(size_t n, double* input, double* output, double W, double norm, info -- information structure about optimization */ int CondatChambollePock2_TV(size_t M, size_t N, double*Y, double lambda, double*X, short alg, int maxit, double* info) { - double sigma, tau, theta, gamma, normalizer, tmp; + double sigma, tau, theta, normalizer, tmp; + double gamma = -1; double *U1 = NULL, *U2 = NULL, *Xt = NULL, *Z = NULL; - int i, j; + size_t i, j; #define FREE \ if(U1) free(U1); \ @@ -695,12 +666,12 @@ int CondatChambollePock2_TV(size_t M, size_t N, double*Y, double lambda, double* // Combiner step: Z = Xt + \theta (Xt - X); for ( i = 0 ; i < M*N ; i++ ) Z[i] = Xt[i] + theta * (Xt[i] - X[i]); - // Update algorithm parameters (only in accelerated Chambolle-Pock) - if ( alg == ALG_CHAMBOLLE_POCK_ACC ) { + // Update algorithm parameters (only in accelerated Chambolle-Pock) + if ( alg == ALG_CHAMBOLLE_POCK_ACC ) { tau *= theta; sigma /= theta; theta = 1. / sqrt(1 + 2 * gamma * tau); // 1/sqrt(1 + 2*gamma*tau) - } + } // Compute stopping tolerance before copying X stop = 0; normalizer = 0; @@ -787,7 +758,7 @@ int CondatChambollePock2_TV(size_t M, size_t N, double*Y, double lambda, double* int Yang2_TV(size_t M, size_t N, double*Y, double lambda, double*X, int maxit, double* info) { double *U1 = NULL, *U2 = NULL, *Z1 = NULL, *Z2 = NULL; double rho; - int i, j; + size_t i, j; Workspace *ws = NULL; #define FREE \ @@ -807,7 +778,7 @@ int Yang2_TV(size_t M, size_t N, double*Y, double lambda, double*X, int maxit, d rho = 10; // Alloc memory - int size = (M > N) ? M : N; + int size = int( (M > N) ? M : N ); U1 = (double*)calloc(M*N,sizeof(double)); U2 = (double*)calloc(M*N,sizeof(double)); Z1 = (double*)malloc(sizeof(double)*M*N); @@ -838,7 +809,7 @@ int Yang2_TV(size_t M, size_t N, double*Y, double lambda, double*X, int maxit, d for ( j = 0 ; j < N ; j++ ) ws->in[j] = -1. / rho * U1[j*M+i] + X[j*M+i]; resetWorkspace(ws); - TV(ws->in, lambda/rho, ws->out, NULL, N, 1, ws); + TV(ws->in, lambda/rho, ws->out, NULL, (int) N, 1, ws); // Recover data for ( j = 0 ; j < N ; j++ ) Z1[j*M+i] = ws->out[j]; @@ -849,7 +820,7 @@ int Yang2_TV(size_t M, size_t N, double*Y, double lambda, double*X, int maxit, d // Copy column data to workspace for ( j = 0 ; j < M ; j++ ) ws->in[j] = -1. / rho * U2[i*M+j] + X[i*M+j]; - TV(ws->in, lambda/rho, ws->out, NULL, M, 1, ws); + TV(ws->in, lambda/rho, ws->out, NULL, (int) M, 1, ws); // Recover data memcpy(Z2+i*M, ws->out, sizeof(double)*M); } @@ -907,7 +878,7 @@ int Yang2_TV(size_t M, size_t N, double*Y, double lambda, double*X, int maxit, d int Kolmogorov2_TV(size_t M, size_t N, double*Y, double lambda, double*X, int maxit, double* info) { double sigma, tau, theta, tmp, normalizer; double *U = NULL, *Xprev = NULL, *row=NULL, *rowout=NULL, *TMP=NULL; - int i, j; + size_t i, j; size_t NM = N*M; #define FREE \ @@ -966,7 +937,7 @@ int Kolmogorov2_TV(size_t M, size_t N, double*Y, double lambda, double*X, int ma // Dual prox for cols for ( j = 0 ; j < NM ; j+=M ) { // Normal prox - TV(TMP+j, lambda/sigma, U+j, NULL, M, 1, NULL); + TV(TMP+j, lambda/sigma, U+j, NULL, (int) M, 1, NULL); // Moreau's identity for ( i = 0 ; i < M ; i++ ) (U+j)[i] = sigma*((TMP+j)[i] - (U+j)[i]); @@ -988,7 +959,7 @@ int Kolmogorov2_TV(size_t M, size_t N, double*Y, double lambda, double*X, int ma for ( j = 0 ; j < N ; j++ ) row[j] = TMP[j*M+i]; // Prox operator - TV(row, lambda/(1.+1./tau), rowout, NULL, N, 1, NULL); + TV(row, lambda/(1.+1./tau), rowout, NULL, (int) N, 1, NULL); // Recover output for ( j = 0 ; j < N ; j++ ) X[j*M+i] = rowout[j]; diff --git a/src/TVL1Wopt.cpp b/src/TVL1Wopt.cpp index afbe566..0f507bf 100644 --- a/src/TVL1Wopt.cpp +++ b/src/TVL1Wopt.cpp @@ -36,10 +36,12 @@ */ int PN_TV1_Weighted(double *y,double *lambda,double *x,double *info,int n,double sigma,Workspace *ws){ int i,ind,nI,recomp,found,iters,nn=n-1; - double lambdaMax,tmp,fval0,fval1,gRd,delta,grad0,stop,stopPrev,improve,rhs,maxStep,prevDelta; + double lambdaMax,tmp,fval0,fval1,gRd,delta,stop,stopPrev,improve,rhs,prevDelta; + double grad0 = 0.; + double maxStep = -DBL_MAX; double *w=NULL,*g=NULL,*d=NULL,*aux=NULL,*aux2=NULL; int *inactive=NULL; - lapack_int one=1,rc,nnp=nn,nIp; + lapack_int nnp=nn,nIp; /* Macros */ #define GRAD2GAP(g,w,gap,i) \ @@ -111,9 +113,15 @@ int PN_TV1_Weighted(double *y,double *lambda,double *x,double *info,int n,double aux2[i] = -1; } aux[nn-1] = 2; +#ifdef PROXTV_USE_LAPACK + lapack_int one = 1; + lapack_int rc; dpttrf_(&nnp,aux,aux2,&rc); /* Solve Choleski-like linear system to obtain unconstrained solution */ dpttrs_(&nnp, &one, aux, aux2, w, &nnp, &rc); +#else + dpttrf_plus_dpttrs_eigen(&nnp, aux, aux2, w); +#endif /* above assume we solved DD'u = Dy */ /* we wanted to solve DD'Wu = Dy; so now obtain u by dividing by W */ @@ -191,17 +199,18 @@ int PN_TV1_Weighted(double *y,double *lambda,double *x,double *info,int n,double #endif /* Factorize reduced Hessian */ nIp = nI; - dpttrf_(&nIp,aux,aux2,&rc); - #ifdef DEBUG - fprintf(DEBUG_FILE,"c=["); for(i=0;ifirst = pb->segments; \ pb->last = pb->first-1; + /******************************************** Taut-string functions ********************************************/ @@ -289,16 +290,18 @@ int classicTautString_TV1_offset(double *signal, int n, double lam, double *prox // Iterate along the signal length Segment *saux; - int i, iaux; + size_t i, iaux; double *pwriter = prox; for ( i = 1 ; i < n-1 ; i++ ) { // Update majorant segment.incx = 1; segment.slope = segment.incy = signal[i]; +#pragma warning(suppress: 4244) concavemajorantadd(majorant, dirsegment, saux, iaux); // Update minorant segment.incx = 1; segment.slope = segment.incy = signal[i]; +#pragma warning(suppress: 4244) convexminorantadd(minorant, dirsegment, saux, iaux); // Update last explored point lastexplored.x++; @@ -316,10 +319,12 @@ int classicTautString_TV1_offset(double *signal, int n, double lam, double *prox // Update majorant with last segment segment.incx = 1; segment.slope = segment.incy = signal[n-1] + lam; +#pragma warning(suppress: 4244) concavemajorantadd(majorant, dirsegment, saux, iaux); // Update minorant with last segment segment.incx = 1; segment.slope = segment.incy = signal[n-1] - lam; +#pragma warning(suppress: 4244) convexminorantadd(minorant, dirsegment, saux, iaux); // At this point, because the endpoint of the tube is the same diff --git a/src/TVL2opt.cpp b/src/TVL2opt.cpp index f628c4c..af36ee3 100644 --- a/src/TVL2opt.cpp +++ b/src/TVL2opt.cpp @@ -36,7 +36,7 @@ int more_TV2(double *y,double lambda,double *x,double *info,int n){ int nn=n-1,i; double stop,tmp,lam,pNorm,qNorm,pNormSq,dist; double *Dy,*alpha,*beta,*minus,*p,*aux; - lapack_int one=1,rc,nnp=nn; + lapack_int nnp=nn; /* Macros */ @@ -101,13 +101,19 @@ int more_TV2(double *y,double lambda,double *x,double *info,int n){ for(i=0;i 1 if(lambdaStep < LAMBDA_STEPS_TVLP){ lambdaCurrent = pow(10, log10(lambdaIni) + lambdaStep * log10(LAMBDA_REDUCTION_TVLP) / ((double)(LAMBDA_STEPS_TVLP-1)) ); GRAD2GAP(w,g,stop,lambdaCurrent,p,n,i,tmp) stuck = 0; bestdual = DBL_MAX; } else break; + #else + break; + #endif } #ifdef TIMING @@ -297,7 +303,7 @@ int OGP_TVp(double *y,double lambda,double *x,double *info,int n,double p,Worksp double q,tmp,stop,dual,bestdual,lambdaMax,lambdaIni,lambdaCurrent,mu,musqrt,beta; int iter,stuck,nn,i,lambdaStep; Workspace *wsinner=NULL; - lapack_int one=1,rc,nnp; + lapack_int nnp; /* Problem constants */ #define L 4 // Lipschitz constant @@ -389,10 +395,15 @@ int OGP_TVp(double *y,double lambda,double *x,double *info,int n,double p,Worksp } aux[nn-1] = 2; nnp=nn; +#ifdef PROXTV_USE_LAPACK + lapack_int one=1; + lapack_int rc; dpttrf_(&nnp,aux,aux2,&rc); /* Solve Choleski-like linear system to obtain unconstrained solution */ dpttrs_(&nnp, &one, aux, aux2, w, &nnp, &rc); - +#else + dpttrf_plus_dpttrs_eigen(&nnp, aux, aux2, w); +#endif /* Compute maximum effective penalty */ lambdaMax = LPnorm(w, nn, q); @@ -514,12 +525,16 @@ int OGP_TVp(double *y,double lambda,double *x,double *info,int n,double p,Worksp while(stop < STOP_TVLP){ /* If met, move to following lambda step, if any */ lambdaStep++; + #if LAMBDA_STEPS_TVLP > 1 if(lambdaStep < LAMBDA_STEPS_TVLP){ lambdaCurrent = pow(10, log10(lambdaIni) + lambdaStep * log10(LAMBDA_REDUCTION_TVLP) / ((double)(LAMBDA_STEPS_TVLP-1)) ); GRAD2GAP(w,g,stop,lambdaCurrent,p,n,i,tmp) stuck = 0; bestdual = DBL_MAX; } else break; + #else + break; + #endif } #ifdef TIMING @@ -582,10 +597,10 @@ int OGP_TVp(double *y,double lambda,double *x,double *info,int n,double p,Worksp */ int FISTA_TVp(double *y,double lambda,double *x,double *info,int n,double p,Workspace *ws){ double *w=NULL,*aux=NULL,*aux2=NULL,*z=NULL,*wpre=NULL,*g; - double q,tmp,stop,dual,bestdual,lambdaMax,lambdaIni,lambdaCurrent,mu,musqrt,beta; + double q,tmp,stop,dual,bestdual,lambdaMax,lambdaIni,lambdaCurrent,beta; int iter,stuck,nn,i,lambdaStep; Workspace *wsinner=NULL; - lapack_int one=1,rc,nnp; + lapack_int nnp; /* Problem constants */ #define L 4 // Lipschitz constant @@ -679,10 +694,15 @@ int FISTA_TVp(double *y,double lambda,double *x,double *info,int n,double p,Work } aux[nn-1] = 2; nnp=nn; +#ifdef PROXTV_USE_LAPACK + lapack_int one=1; + lapack_int rc; dpttrf_(&nnp,aux,aux2,&rc); /* Solve Choleski-like linear system to obtain unconstrained solution */ dpttrs_(&nnp, &one, aux, aux2, w, &nnp, &rc); - +#else + dpttrf_plus_dpttrs_eigen(&nnp, aux, aux2, w); +#endif /* Compute maximum effective penalty */ lambdaMax = LPnorm(w, nn, q); @@ -802,12 +822,16 @@ int FISTA_TVp(double *y,double lambda,double *x,double *info,int n,double p,Work while(stop < STOP_TVLP){ /* If met, move to following lambda step, if any */ lambdaStep++; + #if LAMBDA_STEPS_TVLP > 1 if(lambdaStep < LAMBDA_STEPS_TVLP){ lambdaCurrent = pow(10, log10(lambdaIni) + lambdaStep * log10(LAMBDA_REDUCTION_TVLP) / ((double)(LAMBDA_STEPS_TVLP-1)) ); GRAD2GAP(w,g,stop,lambdaCurrent,p,n,i,tmp) stuck = 0; bestdual = DBL_MAX; } else break; + #else + break; + #endif } #ifdef TIMING @@ -872,7 +896,7 @@ int FW_TVp(double *y,double lambda,double *x,double *info,int n,double p,Workspa int nn, iter, i; double *w=NULL, *aux=NULL, *aux2=NULL, *g=NULL; double stop, q, lambdaMax, gap, step, den, tmp, fval, fvalPrev; - lapack_int one=1,rc,nnp; + lapack_int nnp; /* Gradient to dual gap, considering also the special cases p~=1 and p~=inf */ #define GRAD2GAP(w,g,gap,lambda,p,n,i,tmp) \ @@ -943,9 +967,16 @@ int FW_TVp(double *y,double lambda,double *x,double *info,int n,double p,Workspa } aux[nn-1] = 2; nnp=nn; + +#ifdef PROXTV_USE_LAPACK + lapack_int one = 1; + lapack_int rc; dpttrf_(&nnp,aux,aux2,&rc); /* Solve Choleski-like linear system to obtain unconstrained solution */ dpttrs_(&nnp, &one, aux, aux2, w, &nnp, &rc); +#else + dpttrf_plus_dpttrs_eigen(&nnp, aux, aux2, w); +#endif /* Compute maximum effective penalty */ lambdaMax = LPnorm(w, nn, q); @@ -1111,9 +1142,9 @@ int FW_TVp(double *y,double lambda,double *x,double *info,int n,double p,Workspa int GPFW_TVp(double *y,double lambda,double *x,double *info,int n,double p,Workspace *ws) { double *w=NULL,*aux=NULL,*aux2=NULL,*g=NULL; double q,tmp,stop,dual,bestdual,lambdaMax,num,den,step; - int iter,stuck,nn,i,lambdaStep,cycle; + int iter,stuck,nn,i,cycle; Workspace *wsinner=NULL; - lapack_int one=1,rc,nnp; + lapack_int nnp; /* Problem constants */ #define Linv 0.25 // Inverse of Lipschitz constant @@ -1190,10 +1221,16 @@ int GPFW_TVp(double *y,double lambda,double *x,double *info,int n,double p,Works } aux[nn-1] = 2; nnp=nn; + +#ifdef PROXTV_USE_LAPACK + lapack_int one = 1; + lapack_int rc; dpttrf_(&nnp,aux,aux2,&rc); /* Solve Choleski-like linear system to obtain unconstrained solution */ dpttrs_(&nnp, &one, aux, aux2, w, &nnp, &rc); - +#else + dpttrf_plus_dpttrs_eigen(&nnp, aux, aux2, w); +#endif /* Compute maximum effective penalty */ lambdaMax = LPnorm(w, nn, q); diff --git a/src/TVNDopt.cpp b/src/TVNDopt.cpp index bbde227..bda1c66 100644 --- a/src/TVNDopt.cpp +++ b/src/TVNDopt.cpp @@ -161,7 +161,7 @@ int PD_TV(double *y,double *lambdas,double *norms,double *dims,double *x,double } /* Parallelize */ - #pragma omp parallel shared(ws,nSlices,ns,incs,x,p,lambdas,z,norms,npen,dims,DEBUG_FILE) private(d,i,j,k,idx1,idx2) default(none) + #pragma omp parallel shared(ws,nSlices,ns,incs,x,p,lambdas,z,norms,npen,dims) private(d,i,j,k,idx1,idx2) default(none) { /* Get thread number */ int id = omp_get_thread_num(); @@ -169,10 +169,7 @@ int PD_TV(double *y,double *lambdas,double *norms,double *dims,double *x,double /* Prox step for every penalty term */ for(i=0;iwarm = 0; @@ -187,16 +184,6 @@ int PD_TV(double *y,double *lambdas,double *norms,double *dims,double *x,double for(k=0,idx2=0 ; kin[k] = z[i][idx1+idx2]; - #ifdef DEBUG - { - int dbgi; - fprintf(DEBUG_FILE,"Slice %d: ",j); - for(dbgi=0;dbgiin[dbgi]); - fprintf(DEBUG_FILE,"\n"); - } - #endif - /* Apply 1-dimensional solver */ resetWorkspace(wsi); TV(wsi->in, lambdas[i], wsi->out, NULL, ns[d], norms[i], wsi); @@ -237,9 +224,6 @@ int PD_TV(double *y,double *lambdas,double *norms,double *dims,double *x,double /* Termination check */ if(iters >= MAX_ITERS_PD){ - #ifdef DEBUG - fprintf(DEBUG_FILE,"(PD_TV) WARNING: maximum number of iterations reached (%d).\n",MAX_ITERS_PD); - #endif if(info) info[INFO_RC] = RC_ITERS; } else if(info) info[INFO_RC] = RC_OK; @@ -397,7 +381,7 @@ int PDR_TV(double *y,double *lambdas,double *norms,double *dims,double *x,double } /* Parallelize */ -#pragma omp parallel shared(ws,nSlices,ns,incs,x,p,q,lambdas,z,norms,npen,dims,DEBUG_FILE) private(d,i,j,k,idx1,idx2) default(none) +#pragma omp parallel shared(ws,nSlices,ns,incs,x,p,q,lambdas,z,norms,npen,dims) private(d,i,j,k,idx1,idx2) default(none) { /* Get thread number */ int id = omp_get_thread_num(); @@ -405,10 +389,7 @@ int PDR_TV(double *y,double *lambdas,double *norms,double *dims,double *x,double /* Prox step for every penalty term */ for(i=0;iwarm = 0; @@ -423,16 +404,6 @@ int PDR_TV(double *y,double *lambdas,double *norms,double *dims,double *x,double for(k=0,idx2=0 ; kin[k] = z[i][idx1+idx2]; - #ifdef DEBUG - { - int dbgi; - fprintf(DEBUG_FILE,"Slice %d: ",j); - for(dbgi=0;dbgiin[dbgi]); - fprintf(DEBUG_FILE,"\n"); - } - #endif - /* Apply 1-dimensional solver */ resetWorkspace(wsi); if(norms[i] == 1) @@ -485,9 +456,6 @@ int PDR_TV(double *y,double *lambdas,double *norms,double *dims,double *x,double /* Termination check */ if(iters >= MAX_ITERS_DR){ - #ifdef DEBUG - fprintf(DEBUG_FILE,"(PDR_TV) WARNING: maximum number of iterations reached (%d).\n",MAX_ITERS_DR); - #endif if(info) info[INFO_RC] = RC_ITERS; } else if(info) info[INFO_RC] = RC_OK; @@ -596,7 +564,7 @@ double TVval(double *x,double *lambdas,double *norms,double *dims,int *ns,int nd /* Parallelize calculation of value */ - #pragma omp parallel shared(x,ws,nSlices,ns,incs,lambdas,norms,npen,dims,DEBUG_FILE) private(d,i,j,k,idx1,idx2) default(none) + #pragma omp parallel shared(x,ws,nSlices,ns,incs,lambdas,norms,npen,dims) private(d,i,j,k,idx1,idx2) default(none) { /* Get thread number */ int id = omp_get_thread_num(); @@ -608,10 +576,7 @@ double TVval(double *x,double *lambdas,double *norms,double *dims,int *ns,int nd /* Value for every penalty term */ for(i=0;i N) ? M : N; size = (O > size) ? O : size; - int totalSize = M*N*O; + size_t size_long = (M > N) ? M : N ; size_long = (O > size_long) ? O : size_long; + int size = int(size_long); + size_t totalSize = M*N*O; U1 = (double*)calloc(totalSize,sizeof(double)); U2 = (double*)calloc(totalSize,sizeof(double)); U3 = (double*)calloc(totalSize,sizeof(double)); @@ -738,7 +705,7 @@ int Yang3_TV(size_t M, size_t N, size_t O, double*Y, double lambda, double*X, in ws->in[k] = -1. / rho * U1[idx] + X[idx]; } resetWorkspace(ws); - TV(ws->in, lambda/rho, ws->out, NULL, M, 1, ws); + TV(ws->in, lambda/rho, ws->out, NULL, (int) M, 1, ws); // Recover data memcpy(Z1 + M * ( i + N * j ), ws->out, sizeof(double)*M); } @@ -753,7 +720,7 @@ int Yang3_TV(size_t M, size_t N, size_t O, double*Y, double lambda, double*X, in ws->in[i] = -1. / rho * U2[idx] + X[idx]; } resetWorkspace(ws); - TV(ws->in, lambda/rho, ws->out, NULL, N, 1, ws); + TV(ws->in, lambda/rho, ws->out, NULL, (int) N, 1, ws); // Recover data for ( i = 0 ; i < N ; i++ ) { idx = k + M * ( i + N * j ); @@ -771,7 +738,7 @@ int Yang3_TV(size_t M, size_t N, size_t O, double*Y, double lambda, double*X, in ws->in[j] = -1. / rho * U3[idx] + X[idx]; } resetWorkspace(ws); - TV(ws->in, lambda/rho, ws->out, NULL, O, 1, ws); + TV(ws->in, lambda/rho, ws->out, NULL, (int) O, 1, ws); // Recover data for ( j = 0 ; j < O ; j++ ) { idx = k + M * ( i + N * j ); diff --git a/src/condat_fast_tv.cpp b/src/condat_fast_tv.cpp index 8b17065..ed0c25d 100644 --- a/src/condat_fast_tv.cpp +++ b/src/condat_fast_tv.cpp @@ -195,7 +195,7 @@ void TV1D_denoise_tautstring(double* input, double* output, int width, const dou z[c+i]=y_low[index[c+i]=index_low[s_low+i]]; c = c + c_low-s_low; int j=0; - float a; + floattype a; i=1; while (i<=c) { a = (z[i]-z[i-1]) / (index[i]-index[i-1]); diff --git a/src/general.h b/src/general.h index 96d7d98..bd3efd8 100644 --- a/src/general.h +++ b/src/general.h @@ -15,30 +15,17 @@ #include #include -/* Mex and LAPACK includes */ #ifdef NOMATLAB -#undef lapack_int -#define lapack_int int -extern "C" { - void dpttrs_(lapack_int* n, lapack_int* nrhs, const double* d, const double* e, double* b, lapack_int* ldb, - lapack_int *info ); - void dpttrf_( lapack_int* n, double* d, double* e, lapack_int *info ); -} inline double mxGetInf() { return INFINITY; } -#else -#include "mex.h" -#include "lapack.h" -#include "matrix.h" -#define lapack_int ptrdiff_t #endif +#include "lapackFunctionsWrap.h" + /* Uncomment to print debug messages to a debug file */ //#define DEBUG /* Choose here the name of the debug file */ #ifdef DEBUG static FILE* DEBUG_FILE = fopen("debug.tmp","w"); -#else - static FILE* DEBUG_FILE = NULL; #endif #define DEBUG_N 10 /* Maximum vector length to print in debug messages */ diff --git a/src/lapackFunctionsWrap.cpp b/src/lapackFunctionsWrap.cpp new file mode 100644 index 0000000..ec719a0 --- /dev/null +++ b/src/lapackFunctionsWrap.cpp @@ -0,0 +1,39 @@ +/** + Definitions of functions wrapping lapack. + Declarations depends on compilations options PROXTV_USE_LAPACK + + @author Pablo Hernandez-Cerdan +*/ + +#include "lapackFunctionsWrap.h" + +#ifndef PROXTV_USE_LAPACK // USE_EIGEN +#include +void dpttrf_plus_dpttrs_eigen( lapack_int* n, double* d, double* e, double *b) +{ + using EigenMatrix = Eigen::MatrixXd; + using EigenVector = Eigen::VectorXd; + using EigenVectorMap = Eigen::Map; + // Eigen has to create the full matrix from the diagonal d and subdiagonal e + int mSize = *n; + EigenMatrix eigenM(mSize,mSize); + EigenVectorMap diag(d, mSize); + EigenVectorMap subAndUpperDiag(e, mSize - 1); + EigenVectorMap inputB_outputX(b, mSize); + // Populate matrix + eigenM.diagonal() = diag; + eigenM.diagonal( 1) = subAndUpperDiag; + eigenM.diagonal(-1) = subAndUpperDiag; + + // Factorize using ldlt (ldl is also possible, faster but less accurate) + // A = LDL' + Eigen::LDLT factorization(eigenM); + // A*X = b + // This modifies the input/output pointer: b + inputB_outputX = factorization.solve(inputB_outputX); + // This modifies the input/output pointers: d and e + EigenMatrix factorized = factorization.matrixLDLT(); + diag = factorized.diagonal(); + subAndUpperDiag = factorized.diagonal(-1); +} +#endif diff --git a/src/lapackFunctionsWrap.h b/src/lapackFunctionsWrap.h new file mode 100644 index 0000000..368b048 --- /dev/null +++ b/src/lapackFunctionsWrap.h @@ -0,0 +1,153 @@ +/** + Lapack functions wrapping + This file creates a NAME_wrap layer on top of lapack functions with name: NAME + The wrap function will call lapack or eigen functions depending on the compilation options + WITH LAPACK: + Declare extern NAME_ functions that will be provided by LAPACKE + WITH EIGEN: + Declare NAME_eigen, equivalent to lapack but using EIGEN. + + @author Pablo Hernandez-Cerdan +*/ +#ifdef NOMATLAB +#undef lapack_int +#define lapack_int int +#else // WITH_MATLAB +#include "mex.h" +#ifdef PROXTV_USE_LAPACK +#include "lapack.h" +#endif +#define lapack_int ptrdiff_t +#endif + +#ifdef PROXTV_USE_LAPACK +extern "C" { + /** + * DPTTRS solves a tridiagonal system of the form + * A * X = B + * using the L*D*L' factorization of A computed by DPTTRF. D is a + * diagonal matrix specified in the vector D, L is a unit bidiagonal + * matrix whose subdiagonal is specified in the vector E, and X and B + * are N by NRHS matrices. + * + * Arguments + * ========= + * + * @param n (input) INTEGER + * The order of the tridiagonal matrix A. N >= 0. + * + * @param nrhs (input) INTEGER + * The number of right hand sides, i.e., the number of columns + * of the matrix B. NRHS >= 0. + * + * @param d (input) DOUBLE PRECISION array, dimension (N) + * The n diagonal elements of the diagonal matrix D from the + * L*D*L' factorization of A. + * + * @param e (input) DOUBLE PRECISION array, dimension (N-1) + * The (n-1) subdiagonal elements of the unit bidiagonal factor + * L from the L*D*L' factorization of A. E can also be regarded + * as the superdiagonal of the unit bidiagonal factor U from the + * factorization A = U'*D*U. + * + * @param b (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) + * On entry, the right hand side vectors B for the system of + * linear equations. + * On exit, the solution vectors, X. + * + * @param ldb (input) INTEGER + * The leading dimension of the array B. LDB >= max(1,N). + * + * @param info (output) INTEGER + * = 0: successful exit + * < 0: if INFO = -k, the k-th argument had an illegal value + * + */ + void dpttrs_(lapack_int* n, lapack_int* nrhs, const double* d, + const double* e, double* b, lapack_int* ldb, lapack_int *info ); + /** + * DPTTRF computes the L*D*L' factorization of a real symmetric + * positive definite tridiagonal matrix A. The factorization may also + * be regarded as having the form A = U'*D*U. + * + * Arguments + * ========= + * + * @param n (input) INTEGER + * The order of the matrix A. N >= 0. + * + * @param d (input/output) DOUBLE PRECISION array, dimension (N) + * On entry, the n diagonal elements of the tridiagonal matrix + * A. On exit, the n diagonal elements of the diagonal matrix + * D from the L*D*L' factorization of A. + * + * @param e (input/output) DOUBLE PRECISION array, dimension (N-1) + * On entry, the (n-1) subdiagonal elements of the tridiagonal + * matrix A. On exit, the (n-1) subdiagonal elements of the + * unit bidiagonal factor L from the L*D*L' factorization of A. + * E can also be regarded as the superdiagonal of the unit + * + * @param b (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) + * On entry, the right hand side vectors B for the system of + * linear equations. + * On exit, the solution vectors, X. + * bidiagonal factor U from the U'*D*U factorization of A. + * + * @param info (output) INTEGER + * = 0: successful exit + * < 0: if INFO = -k, the k-th argument had an illegal value + * > 0: if INFO = k, the leading minor of order k is not + * positive definite; if k < N, the factorization could not + * be completed, while if k = N, the factorization was + * completed, but D(N) <= 0. + */ + void dpttrf_( lapack_int* n, double* d, double* e, lapack_int *info ); +} + +#else // USE_EIGEN + +/** + * DPTTRS solves a tridiagonal system of the form + * A * X = B + * using the L*D*L' factorization of A computed by DPTTRF. D is a + * diagonal matrix specified in the vector D, L is a unit bidiagonal + * matrix whose subdiagonal is specified in the vector E, and X and B + * are N by NRHS matrices. + * + * DPTTRF computes the L*D*L' factorization of a real symmetric + * positive definite tridiagonal matrix A. The factorization may also + * be regarded as having the form A = U'*D*U. + * + * Arguments + * ========= + * + * @param n (input) INTEGER + * The order of the matrix A. N >= 0. + * + * @param d (input) DOUBLE PRECISION array, dimension (N) + * On entry, the n diagonal elements of the tridiagonal matrix + * A. On exit, the n diagonal elements of the diagonal matrix + * D from the L*D*L' factorization of A. + * + * @param e (input) DOUBLE PRECISION array, dimension (N-1) + * On entry, the (n-1) subdiagonal elements of the tridiagonal + * matrix A. On exit, the (n-1) subdiagonal elements of the + * unit bidiagonal factor L from the L*D*L' factorization of A. + * E can also be regarded as the superdiagonal of the unit + * + * @param b (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) + * On entry, the right hand side vectors B for the system of + * linear equations. + * On exit, the solution vectors, X. + * + * @param info (output) INTEGER ----- UNUSED ----- + * = 0: successful exit + * < 0: if INFO = -k, the k-th argument had an illegal value + * > 0: if INFO = k, the leading minor of order k is not + * positive definite; if k < N, the factorization could not + * be completed, while if k = N, the factorization was + * completed, but D(N) <= 0. + */ +void dpttrf_plus_dpttrs_eigen( lapack_int* n, double* d, double* e, double* b); + +#endif diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt new file mode 100644 index 0000000..7c02bbf --- /dev/null +++ b/test/CMakeLists.txt @@ -0,0 +1,32 @@ +set(test_names + test_tv1_2d + #other_test + ) + +add_executable(test_tv1_2d "test_tv1_2d.cpp") +add_test(NAME tv1_2d COMMAND test_tv1_2d) + +if(proxTV_USE_LAPACK) + find_package(Eigen3) + if (TARGET Eigen3::Eigen) + list(APPEND test_names test_use_eigen) + add_executable(test_use_eigen "test_use_eigen.cpp") + add_test(NAME use_eigen COMMAND test_use_eigen) + target_link_libraries(test_use_eigen PUBLIC Eigen3::Eigen) + else() + message(WARNING "Eigen3 not found, test comparing lapack and eigen is disabled.") + endif() +endif() +# To add more tests: append a name to test_names, and add_executable/add_test with that name +# Tests can be run with `ctest` in the build directory, +# it accepts regex, and different verbosity: ctest -R tv1 -VV +# add_executable(other_test "a_test_source_file.cpp") +# add_test(NAME a_name_that_makes_sense COMMAND other_test) + +# Link proxTV and include_directory for each executable. +foreach(test_executable ${test_names}) + target_include_directories(${test_executable} PUBLIC + $) + target_compile_definitions(${test_executable} PUBLIC NOMATLAB) + target_link_libraries(${test_executable} PUBLIC proxTV) +endforeach() diff --git a/test/project_using_proxTV/CMakeLists.txt b/test/project_using_proxTV/CMakeLists.txt new file mode 100644 index 0000000..555ff14 --- /dev/null +++ b/test/project_using_proxTV/CMakeLists.txt @@ -0,0 +1,6 @@ +cmake_minimum_required(VERSION 2.8.12) +project(dummy) +find_package(proxTV) +add_executable(dummy dummy.cpp) +# All the requirements of proxTV: compiler definitions, include_directories, linker flags (INTERFACE) are propagated to the target. +target_link_libraries(dummy PUBLIC proxTV::proxTV ) diff --git a/test/project_using_proxTV/README.md b/test/project_using_proxTV/README.md new file mode 100644 index 0000000..3cdac16 --- /dev/null +++ b/test/project_using_proxTV/README.md @@ -0,0 +1,12 @@ +This a minimal project using proxTV as a third-party. +To compile it, it requires that proxTV is installed somewhere in your system. +The installation folder can be modified when configuring proxTV with the option + + -DCMAKE_INSTALL_PREFIX:PATH="/proxTV_install_folder" + +The default is usually the system folder `/usr`. Please read https://github.com/albarji/proxTV/README.md#c-interface for more info. + +Once installed in your system, if it is not in your system path, provide the folder where the `proxTVConfig.cmake` file is located with the option `-proxTV_DIR`: + + mkdir ~/dummy_project ; cd ~/dummy_project + cmake /path/proxTV/test/project_using_proxTV/ -DproxTV_DIR:PATH="/proxTV_install_folder/lib/cmake/proxTV" diff --git a/test/project_using_proxTV/dummy.cpp b/test/project_using_proxTV/dummy.cpp new file mode 100644 index 0000000..876a299 --- /dev/null +++ b/test/project_using_proxTV/dummy.cpp @@ -0,0 +1,26 @@ +#include "TVopt.h" + + +int main() +{ +// int DR2_TV(size_t M, size_t N, double*unary, double W1, double W2, double norm1, double norm2, double*s, int nThreads, int maxit, double* info); + + size_t M = 1; + size_t N = 1; + double* unary = new double; + double W1 = 1.0; + double W2 = 1.0; + double norm1 = 1.0; + double norm2 = 1.0; + double* s = new double; + int nThreads = 1; + int maxit = 1; + double* info = new double; + int r = DR2_TV(M, N, unary, W1, W2, norm1, norm2, s, nThreads, maxit, info); + + delete unary; + delete s; + delete info; + + return r; +} diff --git a/test/test_tv1_2d.cpp b/test/test_tv1_2d.cpp new file mode 100644 index 0000000..876a299 --- /dev/null +++ b/test/test_tv1_2d.cpp @@ -0,0 +1,26 @@ +#include "TVopt.h" + + +int main() +{ +// int DR2_TV(size_t M, size_t N, double*unary, double W1, double W2, double norm1, double norm2, double*s, int nThreads, int maxit, double* info); + + size_t M = 1; + size_t N = 1; + double* unary = new double; + double W1 = 1.0; + double W2 = 1.0; + double norm1 = 1.0; + double norm2 = 1.0; + double* s = new double; + int nThreads = 1; + int maxit = 1; + double* info = new double; + int r = DR2_TV(M, N, unary, W1, W2, norm1, norm2, s, nThreads, maxit, info); + + delete unary; + delete s; + delete info; + + return r; +} diff --git a/test/test_use_eigen.cpp b/test/test_use_eigen.cpp new file mode 100644 index 0000000..3acaf3e --- /dev/null +++ b/test/test_use_eigen.cpp @@ -0,0 +1,178 @@ +#include +#include +#include "general.h" + +using EigenMatrix = Eigen::MatrixXd; +using EigenVector = Eigen::VectorXd; + +EigenMatrix oneD() +{ + EigenMatrix m(1,4); + m(0,0) = 3; + m(0,1) = 2.5; + m(0,2) = -1; + m(0,3) = 1.5; + return m; +} + +EigenMatrix lapack_oneD() +{ + auto eigenM = oneD(); + int n = eigenM.cols(); + std::cout << "Cols: " << n << std::endl; + int nn = n ; + int rc = nn; + int nnp = nn; + int one = 1; + double *w = nullptr; + double *aux = nullptr; + double *aux2 = nullptr; + w = (double*)malloc(sizeof(double)*nn); + aux = (double*)malloc(sizeof(double)*nn); + aux2 = (double*)malloc(sizeof(double)*nn); + for(int i=0;i ldltOfM(eigenM); + EigenVector v = ldltOfM.solve(b); + std::cout << v << std::endl; + // EigenVector v = eigenM.colPivHouseholderQr().solve(b); + EigenMatrix ldltMatrix = ldltOfM.matrixLDLT(); + std::cout << ldltMatrix << std::endl; + EigenMatrix lMatrix = ldltOfM.matrixL(); + std::cout << lMatrix << std::endl; + std::cout << "Diagonal of LDLT:" << std::endl; + std::cout << ldltMatrix.diagonal() << std::endl; + std::cout << "SubDiagonal of LDLT:" << std::endl; + std::cout << ldltMatrix.diagonal(-1) << std::endl; + + return v; +} + +/// The idea is that input and output are raw pointers to couple with existing code +EigenMatrix eigen_interface_raw_pointers_oneD() +{ + // Create the typemaps + using EigenMatrixMap = Eigen::Map; + using EigenMatrixReadOnlyMap = Eigen::Map; + using EigenVectorMap = Eigen::Map; + using EigenVectorReadOnlyMap = Eigen::Map; + // Populate values using Eigen (just because easier) + EigenMatrix eigenM(4,4); + { + EigenVector diag(eigenM.cols()); + diag.setConstant(2); + EigenVector subUpperDiag(eigenM.cols()); + subUpperDiag.setConstant(-1); + eigenM.diagonal() = diag; + eigenM.diagonal(1) = subUpperDiag; + eigenM.diagonal(-1) = subUpperDiag; + } + std::cout << "Matrix eigenM:\n" << eigenM << std::endl; + EigenVector b(eigenM.cols()); + EigenMatrix oneM = oneD(); + for (int i = 0; i < oneM.cols(); ++i) { + b[i] = oneM(0, i+1) - oneM(0, i); /* Dy */ + } + // Set input raw pointers + int nn = eigenM.cols(); + double *w = (double*)malloc(sizeof(double)*nn); + for(int i=0;i