diff --git a/src/external/CoMISo/CHANGELOG b/src/external/CoMISo/CHANGELOG
new file mode 100644
index 000000000..1b46e752e
--- /dev/null
+++ b/src/external/CoMISo/CHANGELOG
@@ -0,0 +1,38 @@
+#============================================================================================
+
+CoMISo 1.0-rc1:
+- Initial release
+
+CoMISo 1.1:
+- Only 2 external dependencies: header libraries GMM++ and Eigen3
+- Efficient re-solve with updated (system and/or constraints) right hand sides
+
+
+Features:
+ - /NSolver -- simple non-linear solver c++ interfaces:
+ * Gurobi [1],
+ * IPOPT [2],
+ * TAO/PETSc [3],
+ * CPLEX [4],
+ * an the inhouse Newton Solver and more!
+
+ - /EigenSolver -- interface for large-scale eigenvalue problems ARPACK [5]
+
+ - /Solver -- several interfaces for common linear system solvers:
+ * sparse Cholesky (Cholmod, LDLT)
+ * SparseQR [6]
+ * TAUCS [7]
+ * UMFPACK [8]
+ and of course the Constrained and Constrained Mixed Integer Solvers (CoMISo)
+
+
+
+References:
+[1] www.gurobi.com
+[2] https://projects.coin-or.org/Ipopt
+[3] http://www.mcs.anl.gov/research/projects/tao/
+[4] http://www-01.ibm.com/software/integration/optimization/cplex-optimizer/
+[5] http://www.caam.rice.edu/software/ARPACK/
+[6] http://www.cise.ufl.edu/research/sparse/SPQR/
+[7] http://www.tau.ac.il/~stoledo/taucs/
+[8] http://www.cise.ufl.edu/research/sparse/umfpack/
diff --git a/src/external/CoMISo/CMakeLists.txt b/src/external/CoMISo/CMakeLists.txt
new file mode 100644
index 000000000..131d6bb78
--- /dev/null
+++ b/src/external/CoMISo/CMakeLists.txt
@@ -0,0 +1,371 @@
+cmake_minimum_required (VERSION 2.6)
+
+project(CoMISo)
+
+# add our macro directory to cmake search path
+set (CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}/cmake)
+
+include (ACGCommon)
+
+acg_qt4 ()
+# change to 0 if QT should not be used
+set( WANT_COMISO_QT 1 )
+if( QT4_FOUND)
+ #message( WARNING " QT4 FOUND" )
+ if( WANT_COMISO_QT )
+ add_definitions (-DQT4_FOUND)
+ # message( WARNING " USING QT4" )
+ endif ()
+ set (COMISO_QT4_CONFIG_FILE_SETTINGS "#define COMISO_QT4_AVAILABLE 1" )
+else()
+ set (COMISO_QT4_CONFIG_FILE_SETTINGS "#define COMISO_QT4_AVAILABLE 0" )
+endif ()
+
+acg_get_version ()
+
+include (ACGOutput)
+
+set(COMISO_INCLUDE_DIRECTORIES "")
+set(COMISO_LINK_DIRECTORIES "")
+set(COMISO_LINK_LIBRARIES "")
+
+FIND_PACKAGE( Boost 1.42.0 COMPONENTS system filesystem regex )
+if(Boost_FOUND)
+ set (COMISO_BOOST_CONFIG_FILE_SETTINGS "#define COMISO_BOOST_AVAILABLE 1" )
+ list( APPEND COMISO_INCLUDE_DIRECTORIES ${Boost_INCLUDE_DIRS} )
+# list( APPEND COMISO_LINK_DIRECTORIES ${Boost_LIBRARY_DIR} )
+ list( APPEND COMISO_LINK_LIBRARIES ${Boost_LIBRARIES} )
+else()
+ set (COMISO_BOOST_CONFIG_FILE_SETTINGS "#define COMISO_BOOST_AVAILABLE 0" )
+ message (STATUS "Boost not found!")
+endif ()
+
+
+find_package (GMM)
+if (GMM_FOUND)
+ set (COMISO_GMM_CONFIG_FILE_SETTINGS "#define COMISO_GMM_AVAILABLE 1" )
+ list( APPEND COMISO_INCLUDE_DIRECTORIES ${GMM_INCLUDE_DIR} )
+else()
+ set (COMISO_GMM_CONFIG_FILE_SETTINGS "#define COMISO_GMM_AVAILABLE 0" )
+ message (FATAL_ERROR "GMM not found!")
+endif ()
+
+# We require cgal with its blas on windows
+ find_package(CGAL)
+if (CGAL_FOUND)
+ set (COMISO_CGAL_CONFIG_FILE_SETTINGS "#define COMISO_CGAL_AVAILABLE 1" )
+ list( APPEND COMISO_INCLUDE_DIRECTORIES ${CGAL_INCLUDE_DIR} )
+ list( APPEND COMISO_LINK_DIRECTORIES ${CGAL_LIBRARY_DIR} )
+ list( APPEND COMISO_LINK_LIBRARIES ${CGAL_LIBRARIES} )
+else()
+ set (COMISO_CGAL_CONFIG_FILE_SETTINGS "#define COMISO_CGAL_AVAILABLE 0" )
+ message (STATUS "CGAL not found!")
+endif()
+
+
+find_package (BLAS)
+if (BLAS_FOUND )
+ set (COMISO_BLAS_CONFIG_FILE_SETTINGS "#define COMISO_BLAS_AVAILABLE 1" )
+
+ list( APPEND COMISO_INCLUDE_DIRECTORIES ${BLAS_INCLUDE_DIRS} )
+ list( APPEND COMISO_LINK_DIRECTORIES ${BLAS_LIBRARY_DIRS} )
+ list( APPEND COMISO_LINK_LIBRARIES ${BLAS_LIBRARIES} )
+else()
+ set (COMISO_BLAS_CONFIG_FILE_SETTINGS "#define COMISO_BLAS_AVAILABLE 0" )
+ message (STATUS "BLAS not found!")
+endif ()
+
+find_package (ADOLC)
+if (ADOLC_FOUND)
+ set (COMISO_ADOLC_CONFIG_FILE_SETTINGS "#define COMISO_ADOLC_AVAILABLE 1" )
+ list( APPEND COMISO_INCLUDE_DIRECTORIES ${ADOLC_INCLUDE_DIR} )
+ list( APPEND COMISO_LINK_DIRECTORIES ${ADOLC_LIBRARY_DIR} )
+ list( APPEND COMISO_LINK_LIBRARIES ${ADOLC_LIBRARIES} )
+else ()
+ set (COMISO_ADOLC_CONFIG_FILE_SETTINGS "#define COMISO_ADOLC_AVAILABLE 0" )
+ message (STATUS "ADOLC not found!")
+endif ()
+
+find_package (SUITESPARSE)
+if (SUITESPARSE_FOUND )
+ set (COMISO_SUITESPARSE_CONFIG_FILE_SETTINGS "#define COMISO_SUITESPARSE_AVAILABLE 1" )
+
+ list( APPEND COMISO_INCLUDE_DIRECTORIES ${SUITESPARSE_INCLUDE_DIRS} )
+ list( APPEND COMISO_LINK_DIRECTORIES ${SUITESPARSE_LIBRARY_DIRS} )
+ list( APPEND COMISO_LINK_LIBRARIES ${SUITESPARSE_LIBRARIES} )
+else ()
+ message (STATUS "SUITESPARSE not found!")
+ set (COMISO_SUITESPARSE_CONFIG_FILE_SETTINGS "#define COMISO_SUITESPARSE_AVAILABLE 0" )
+endif ()
+# special handling, since spqr is incorrect in several distributions
+if(SUITESPARSE_SPQR_VALID)
+ set (COMISO_SUITESPARSE_SPQR_CONFIG_FILE_SETTINGS "#define COMISO_SUITESPARSE_SPQR_AVAILABLE 1" )
+else()
+ message (STATUS "SUITESPARSE SPQR seems to be invalid!")
+ set (COMISO_SUITESPARSE_SPQR_CONFIG_FILE_SETTINGS "#define COMISO_SUITESPARSE_SPQR_AVAILABLE 0" )
+endif()
+
+find_package (MPI)
+if (MPI_FOUND )
+ set (COMISO_MPI_CONFIG_FILE_SETTINGS "#define COMISO_MPI_AVAILABLE 1" )
+ list( APPEND COMISO_INCLUDE_DIRECTORIES ${MPI_INCLUDE_PATH} )
+ list( APPEND COMISO_LINK_LIBRARIES ${MPI_CXX_LIBRARIES} )
+else ()
+ message (STATUS "MPI not found!")
+ set (COMISO_MPI_CONFIG_FILE_SETTINGS "#define COMISO_MPI_AVAILABLE 0" )
+endif ()
+
+find_package (PETSC)
+if (PETSC_FOUND AND MPI_FOUND)
+ set (COMISO_PETSC_CONFIG_FILE_SETTINGS "#define COMISO_PETSC_AVAILABLE 1" )
+ list( APPEND COMISO_INCLUDE_DIRECTORIES ${PETSC_INCLUDE_DIRS} )
+ list( APPEND COMISO_LINK_DIRECTORIES ${PETSC_LIBRARY_DIR} )
+ list( APPEND COMISO_LINK_LIBRARIES ${PETSC_LIBRARY} )
+else ()
+ message (STATUS "PETSC or dependency not found!")
+ set (COMISO_PETSC_CONFIG_FILE_SETTINGS "#define COMISO_PETSC_AVAILABLE 0" )
+endif ()
+
+
+find_package (TAO)
+if (TAO_FOUND AND PETSC_FOUND AND MPI_FOUND)
+ set (COMISO_TAO_CONFIG_FILE_SETTINGS "#define COMISO_TAO_AVAILABLE 1" )
+ list( APPEND COMISO_INCLUDE_DIRECTORIES ${TAO_INCLUDE_DIRS} )
+ list( APPEND COMISO_LINK_DIRECTORIES ${TAO_LIBRARY_DIR} )
+ list( APPEND COMISO_LINK_LIBRARIES ${TAO_LIBRARY} )
+else ()
+ message (STATUS "TAO or dependency not found!")
+ set (COMISO_TAO_CONFIG_FILE_SETTINGS "#define COMISO_TAO_AVAILABLE 0" )
+endif ()
+
+find_package (METIS)
+if (METIS_FOUND )
+ set (COMISO_METIS_CONFIG_FILE_SETTINGS "#define COMISO_METIS_AVAILABLE 1" )
+
+ list( APPEND COMISO_INCLUDE_DIRECTORIES ${METIS_INCLUDE_DIRS} )
+ list( APPEND COMISO_LINK_DIRECTORIES ${METIS_LIBRARY_DIRS} )
+ list( APPEND COMISO_LINK_LIBRARIES ${METIS_LIBRARIES} )
+else()
+ set (COMISO_METIS_CONFIG_FILE_SETTINGS "#define COMISO_METIS_AVAILABLE 0" )
+ message (STATUS "METIS not found!")
+endif ()
+
+find_package (MUMPS)
+if (MUMPS_FOUND )
+ set (COMISO_MUMPS_CONFIG_FILE_SETTINGS "#define COMISO_MUMPS_AVAILABLE 1" )
+ list( APPEND COMISO_INCLUDE_DIRECTORIES ${MUMPS_INCLUDE_DIR} )
+ list( APPEND COMISO_LINK_LIBRARIES ${MUMPS_LIBRARY} )
+else ()
+ message (STATUS "MUMPS not found!")
+ set (COMISO_MUMPS_CONFIG_FILE_SETTINGS "#define COMISO_MUMPS_AVAILABLE 0" )
+endif ()
+
+find_package (IPOPT)
+if (IPOPT_FOUND)
+ set (COMISO_IPOPT_CONFIG_FILE_SETTINGS "#define COMISO_IPOPT_AVAILABLE 1" )
+ list( APPEND COMISO_INCLUDE_DIRECTORIES ${IPOPT_INCLUDE_DIR} )
+ list( APPEND COMISO_LINK_DIRECTORIES ${IPOPT_LIBRARY_DIR} )
+ list( APPEND COMISO_LINK_LIBRARIES ${IPOPT_LIBRARY} )
+ if ( IPOPT_HSL_LIBRARY_DIR )
+ set (COMISO_HSL_CONFIG_FILE_SETTINGS "#define COMISO_HSL_AVAILABLE 1" )
+ else ()
+ set (COMISO_HSL_CONFIG_FILE_SETTINGS "#define COMISO_HSL_AVAILABLE 0" )
+ endif()
+else ()
+ message (STATUS "IPOPT or dependency not found!")
+ set (COMISO_IPOPT_CONFIG_FILE_SETTINGS "#define COMISO_IPOPT_AVAILABLE 0" )
+ set (COMISO_HSL_CONFIG_FILE_SETTINGS "#define COMISO_HSL_AVAILABLE 0" )
+endif ()
+
+find_package (EIGEN3)
+if (EIGEN3_FOUND )
+ set (COMISO_EIGEN3_CONFIG_FILE_SETTINGS "#define COMISO_EIGEN3_AVAILABLE 1" )
+ list( APPEND COMISO_INCLUDE_DIRECTORIES ${EIGEN3_INCLUDE_DIR} )
+else ()
+ message (STATUS "EIGEN3 not found!")
+ set (COMISO_EIGEN3_CONFIG_FILE_SETTINGS "#define COMISO_EIGEN3_AVAILABLE 0" )
+endif ()
+
+find_package (Taucs)
+if (TAUCS_FOUND )
+ set (COMISO_TAUCS_CONFIG_FILE_SETTINGS "#define COMISO_TAUCS_AVAILABLE 1" )
+ list( APPEND COMISO_INCLUDE_DIRECTORIES ${TAUCS_INCLUDE_DIR} )
+ list( APPEND COMISO_INCLUDE_DIRECTORIES ${LAPACK_INCLUDE_DIR} )
+ list( APPEND COMISO_LINK_DIRECTORIES ${LAPACK_LIBRARY_DIR} )
+ list( APPEND COMISO_LINK_LIBRARIES ${TAUCS_LIBRARY} )
+ list( APPEND COMISO_LINK_LIBRARIES ${LAPACK_LIBRARIES} )
+else ()
+ message (STATUS "TAUCS not found!")
+ set (COMISO_TAUCS_CONFIG_FILE_SETTINGS "#define COMISO_TAUCS_AVAILABLE 0" )
+endif ()
+
+find_package (GUROBI)
+if (GUROBI_FOUND )
+ set (COMISO_GUROBI_CONFIG_FILE_SETTINGS "#define COMISO_GUROBI_AVAILABLE 1" )
+ list( APPEND COMISO_INCLUDE_DIRECTORIES ${GUROBI_INCLUDE_DIRS} )
+# list( APPEND COMISO_LINK_DIRECTORIES ${GUROBI_LIBRARY_DIR} )
+ list( APPEND COMISO_LINK_LIBRARIES ${GUROBI_LIBRARIES} )
+else ()
+ message (STATUS "GUROBI not found!")
+ set (COMISO_GUROBI_CONFIG_FILE_SETTINGS "#define COMISO_GUROBI_AVAILABLE 0" )
+endif ()
+
+find_package (ARPACK)
+if (ARPACK_FOUND )
+ set (COMISO_ARPACK_CONFIG_FILE_SETTINGS "#define COMISO_ARPACK_AVAILABLE 1" )
+ list( APPEND COMISO_INCLUDE_DIRECTORIES ${ARPACK_INCLUDE_DIR} )
+# list( APPEND COMISO_LINK_DIRECTORIES ${ARPACK_LIBRARY_DIR} )
+ list( APPEND COMISO_LINK_LIBRARIES ${ARPACK_LIBRARY} )
+else ()
+ message (STATUS "ARPACK not found!")
+ set (COMISO_ARPACK_CONFIG_FILE_SETTINGS "#define COMISO_ARPACK_AVAILABLE 0" )
+endif ()
+
+find_package (CPLEX)
+if (CPLEX_FOUND )
+ set (COMISO_CPLEX_CONFIG_FILE_SETTINGS "#define COMISO_CPLEX_AVAILABLE 1" )
+ list( APPEND COMISO_INCLUDE_DIRECTORIES ${CPLEX_INCLUDE_DIRS} )
+ list( APPEND COMISO_LINK_LIBRARIES ${CPLEX_LIBRARIES} )
+
+ #enable c++ support
+ add_definitions(-DIL_STD)
+else ()
+ message (STATUS "CPLEX not found!")
+ set (COMISO_CPLEX_CONFIG_FILE_SETTINGS "#define COMISO_CPLEX_AVAILABLE 0" )
+endif ()
+
+include_directories (
+ ..
+ ${CMAKE_CURRENT_SOURCE_DIR}
+ ${CMAKE_CURRENT_SOURCE_DIR}/../
+ ${CMAKE_CURRENT_BINARY_DIR}
+ ${COMISO_INCLUDE_DIRECTORIES}
+)
+
+# generate dllexport macros on windows
+if (WIN32)
+ add_definitions(-DCOMISODLL)
+ add_definitions(-D_SCL_SECURE_NO_DEPRECATE)
+endif ()
+
+
+link_directories (
+ ${COMISO_LINK_DIRECTORIES}
+)
+
+# source code directories
+set (directories
+ .
+ Solver
+ NSolver
+ EigenSolver
+ Config
+ Utils
+ QtWidgets
+)
+
+if (WIN32)
+ add_definitions(
+ -D_USE_MATH_DEFINES -DNOMINMAX
+ )
+endif ()
+
+# collect all header,source and ui files
+acg_append_files (headers "*.hh" ${directories})
+acg_append_files (sources "*.cc" ${directories})
+acg_append_files (ui "*.ui" ${directories})
+
+
+macro (of_list_filter _list)
+ if (WIN32)
+ foreach (_element ${${_list}})
+ if (_element MATCHES "gnuplot_i\\.(cc|hh)$")
+ list (REMOVE_ITEM ${_list} ${_element})
+ endif ()
+ endforeach ()
+ endif ()
+endmacro ()
+
+of_list_filter ( headers )
+of_list_filter ( sources )
+
+
+# remove template cc files from source file list
+acg_drop_templates (sources)
+
+if( QT4_FOUND)
+# genereate uic and moc targets
+acg_qt4_autouic (uic_targets ${ui})
+acg_qt4_automoc (moc_targets ${headers})
+endif()
+acg_add_library (CoMISo SHARED ${uic_targets} ${sources} ${headers} ${moc_targets})
+
+if (NOT APPLE)
+ target_link_libraries (CoMISo
+ ${QT_LIBRARIES}
+ ${COMISO_LINK_LIBRARIES}
+ )
+else(NOT APPLE)
+ target_link_libraries (CoMISo
+ ${QT_LIBRARIES}
+ ${COMISO_LINK_LIBRARIES}
+ )
+endif(NOT APPLE)
+
+# display results
+acg_print_configure_header (COMISO "CoMISo")
+
+# write config file
+configure_file ("${CMAKE_CURRENT_SOURCE_DIR}/Config/config.hh.in"
+ "${CMAKE_CURRENT_SOURCE_DIR}/Config/config.hh" @ONLY IMMEDIATE)
+
+
+
+#######################################################################
+# Configure the examples last to be sure, that all configure files
+# of the library are already there
+#######################################################################
+
+if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/factored_solver/CMakeLists.txt" )
+ add_subdirectory (Examples/factored_solver)
+endif()
+
+if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/quadratic_solver/CMakeLists.txt" )
+ add_subdirectory (Examples/quadratic_solver)
+endif()
+if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/test2/CMakeLists.txt" )
+ add_subdirectory (Examples/test2)
+endif()
+if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_quadratic_example/CMakeLists.txt" )
+ add_subdirectory (Examples/small_quadratic_example)
+endif()
+if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_factored_example/CMakeLists.txt" )
+ add_subdirectory (Examples/small_factored_example)
+endif()
+if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/super_sparse_matrix/CMakeLists.txt" )
+ add_subdirectory (Examples/super_sparse_matrix)
+endif()
+if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/eigen_solver/CMakeLists.txt" )
+ add_subdirectory (Examples/eigen_solver)
+endif()
+if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_nsolver/CMakeLists.txt" )
+ add_subdirectory (Examples/small_nsolver)
+endif()
+if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_eigenproblem/CMakeLists.txt" )
+ add_subdirectory (Examples/small_eigenproblem)
+endif()
+if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_miqp/CMakeLists.txt" )
+ add_subdirectory (Examples/small_miqp)
+endif()
+if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_nleast_squares/CMakeLists.txt" )
+ add_subdirectory (Examples/small_nleast_squares)
+endif()
+if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_sparseqr/CMakeLists.txt" )
+ add_subdirectory (Examples/small_sparseqr)
+endif()
+if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_quadratic_resolve_example/CMakeLists.txt" )
+ add_subdirectory (Examples/small_quadratic_resolve_example)
+endif()
+if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_cplex_soc/CMakeLists.txt" )
+ add_subdirectory (Examples/small_cplex_soc)
+endif()
diff --git a/src/external/CoMISo/COPYING b/src/external/CoMISo/COPYING
new file mode 100644
index 000000000..10926e87f
--- /dev/null
+++ b/src/external/CoMISo/COPYING
@@ -0,0 +1,675 @@
+ GNU GENERAL PUBLIC LICENSE
+ Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc.
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+ Preamble
+
+ The GNU General Public License is a free, copyleft license for
+software and other kinds of works.
+
+ The licenses for most software and other practical works are designed
+to take away your freedom to share and change the works. By contrast,
+the GNU General Public License is intended to guarantee your freedom to
+share and change all versions of a program--to make sure it remains free
+software for all its users. We, the Free Software Foundation, use the
+GNU General Public License for most of our software; it applies also to
+any other work released this way by its authors. You can apply it to
+your programs, too.
+
+ When we speak of free software, we are referring to freedom, not
+price. Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+them if you wish), that you receive source code or can get it if you
+want it, that you can change the software or use pieces of it in new
+free programs, and that you know you can do these things.
+
+ To protect your rights, we need to prevent others from denying you
+these rights or asking you to surrender the rights. Therefore, you have
+certain responsibilities if you distribute copies of the software, or if
+you modify it: responsibilities to respect the freedom of others.
+
+ For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must pass on to the recipients the same
+freedoms that you received. You must make sure that they, too, receive
+or can get the source code. And you must show them these terms so they
+know their rights.
+
+ Developers that use the GNU GPL protect your rights with two steps:
+(1) assert copyright on the software, and (2) offer you this License
+giving you legal permission to copy, distribute and/or modify it.
+
+ For the developers' and authors' protection, the GPL clearly explains
+that there is no warranty for this free software. For both users' and
+authors' sake, the GPL requires that modified versions be marked as
+changed, so that their problems will not be attributed erroneously to
+authors of previous versions.
+
+ Some devices are designed to deny users access to install or run
+modified versions of the software inside them, although the manufacturer
+can do so. This is fundamentally incompatible with the aim of
+protecting users' freedom to change the software. The systematic
+pattern of such abuse occurs in the area of products for individuals to
+use, which is precisely where it is most unacceptable. Therefore, we
+have designed this version of the GPL to prohibit the practice for those
+products. If such problems arise substantially in other domains, we
+stand ready to extend this provision to those domains in future versions
+of the GPL, as needed to protect the freedom of users.
+
+ Finally, every program is threatened constantly by software patents.
+States should not allow patents to restrict development and use of
+software on general-purpose computers, but in those that do, we wish to
+avoid the special danger that patents applied to a free program could
+make it effectively proprietary. To prevent this, the GPL assures that
+patents cannot be used to render the program non-free.
+
+ The precise terms and conditions for copying, distribution and
+modification follow.
+
+ TERMS AND CONDITIONS
+
+ 0. Definitions.
+
+ "This License" refers to version 3 of the GNU General Public License.
+
+ "Copyright" also means copyright-like laws that apply to other kinds of
+works, such as semiconductor masks.
+
+ "The Program" refers to any copyrightable work licensed under this
+License. Each licensee is addressed as "you". "Licensees" and
+"recipients" may be individuals or organizations.
+
+ To "modify" a work means to copy from or adapt all or part of the work
+in a fashion requiring copyright permission, other than the making of an
+exact copy. The resulting work is called a "modified version" of the
+earlier work or a work "based on" the earlier work.
+
+ A "covered work" means either the unmodified Program or a work based
+on the Program.
+
+ To "propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy. Propagation includes copying,
+distribution (with or without modification), making available to the
+public, and in some countries other activities as well.
+
+ To "convey" a work means any kind of propagation that enables other
+parties to make or receive copies. Mere interaction with a user through
+a computer network, with no transfer of a copy, is not conveying.
+
+ An interactive user interface displays "Appropriate Legal Notices"
+to the extent that it includes a convenient and prominently visible
+feature that (1) displays an appropriate copyright notice, and (2)
+tells the user that there is no warranty for the work (except to the
+extent that warranties are provided), that licensees may convey the
+work under this License, and how to view a copy of this License. If
+the interface presents a list of user commands or options, such as a
+menu, a prominent item in the list meets this criterion.
+
+ 1. Source Code.
+
+ The "source code" for a work means the preferred form of the work
+for making modifications to it. "Object code" means any non-source
+form of a work.
+
+ A "Standard Interface" means an interface that either is an official
+standard defined by a recognized standards body, or, in the case of
+interfaces specified for a particular programming language, one that
+is widely used among developers working in that language.
+
+ The "System Libraries" of an executable work include anything, other
+than the work as a whole, that (a) is included in the normal form of
+packaging a Major Component, but which is not part of that Major
+Component, and (b) serves only to enable use of the work with that
+Major Component, or to implement a Standard Interface for which an
+implementation is available to the public in source code form. A
+"Major Component", in this context, means a major essential component
+(kernel, window system, and so on) of the specific operating system
+(if any) on which the executable work runs, or a compiler used to
+produce the work, or an object code interpreter used to run it.
+
+ The "Corresponding Source" for a work in object code form means all
+the source code needed to generate, install, and (for an executable
+work) run the object code and to modify the work, including scripts to
+control those activities. However, it does not include the work's
+System Libraries, or general-purpose tools or generally available free
+programs which are used unmodified in performing those activities but
+which are not part of the work. For example, Corresponding Source
+includes interface definition files associated with source files for
+the work, and the source code for shared libraries and dynamically
+linked subprograms that the work is specifically designed to require,
+such as by intimate data communication or control flow between those
+subprograms and other parts of the work.
+
+ The Corresponding Source need not include anything that users
+can regenerate automatically from other parts of the Corresponding
+Source.
+
+ The Corresponding Source for a work in source code form is that
+same work.
+
+ 2. Basic Permissions.
+
+ All rights granted under this License are granted for the term of
+copyright on the Program, and are irrevocable provided the stated
+conditions are met. This License explicitly affirms your unlimited
+permission to run the unmodified Program. The output from running a
+covered work is covered by this License only if the output, given its
+content, constitutes a covered work. This License acknowledges your
+rights of fair use or other equivalent, as provided by copyright law.
+
+ You may make, run and propagate covered works that you do not
+convey, without conditions so long as your license otherwise remains
+in force. You may convey covered works to others for the sole purpose
+of having them make modifications exclusively for you, or provide you
+with facilities for running those works, provided that you comply with
+the terms of this License in conveying all material for which you do
+not control copyright. Those thus making or running the covered works
+for you must do so exclusively on your behalf, under your direction
+and control, on terms that prohibit them from making any copies of
+your copyrighted material outside their relationship with you.
+
+ Conveying under any other circumstances is permitted solely under
+the conditions stated below. Sublicensing is not allowed; section 10
+makes it unnecessary.
+
+ 3. Protecting Users' Legal Rights From Anti-Circumvention Law.
+
+ No covered work shall be deemed part of an effective technological
+measure under any applicable law fulfilling obligations under article
+11 of the WIPO copyright treaty adopted on 20 December 1996, or
+similar laws prohibiting or restricting circumvention of such
+measures.
+
+ When you convey a covered work, you waive any legal power to forbid
+circumvention of technological measures to the extent such circumvention
+is effected by exercising rights under this License with respect to
+the covered work, and you disclaim any intention to limit operation or
+modification of the work as a means of enforcing, against the work's
+users, your or third parties' legal rights to forbid circumvention of
+technological measures.
+
+ 4. Conveying Verbatim Copies.
+
+ You may convey verbatim copies of the Program's source code as you
+receive it, in any medium, provided that you conspicuously and
+appropriately publish on each copy an appropriate copyright notice;
+keep intact all notices stating that this License and any
+non-permissive terms added in accord with section 7 apply to the code;
+keep intact all notices of the absence of any warranty; and give all
+recipients a copy of this License along with the Program.
+
+ You may charge any price or no price for each copy that you convey,
+and you may offer support or warranty protection for a fee.
+
+ 5. Conveying Modified Source Versions.
+
+ You may convey a work based on the Program, or the modifications to
+produce it from the Program, in the form of source code under the
+terms of section 4, provided that you also meet all of these conditions:
+
+ a) The work must carry prominent notices stating that you modified
+ it, and giving a relevant date.
+
+ b) The work must carry prominent notices stating that it is
+ released under this License and any conditions added under section
+ 7. This requirement modifies the requirement in section 4 to
+ "keep intact all notices".
+
+ c) You must license the entire work, as a whole, under this
+ License to anyone who comes into possession of a copy. This
+ License will therefore apply, along with any applicable section 7
+ additional terms, to the whole of the work, and all its parts,
+ regardless of how they are packaged. This License gives no
+ permission to license the work in any other way, but it does not
+ invalidate such permission if you have separately received it.
+
+ d) If the work has interactive user interfaces, each must display
+ Appropriate Legal Notices; however, if the Program has interactive
+ interfaces that do not display Appropriate Legal Notices, your
+ work need not make them do so.
+
+ A compilation of a covered work with other separate and independent
+works, which are not by their nature extensions of the covered work,
+and which are not combined with it such as to form a larger program,
+in or on a volume of a storage or distribution medium, is called an
+"aggregate" if the compilation and its resulting copyright are not
+used to limit the access or legal rights of the compilation's users
+beyond what the individual works permit. Inclusion of a covered work
+in an aggregate does not cause this License to apply to the other
+parts of the aggregate.
+
+ 6. Conveying Non-Source Forms.
+
+ You may convey a covered work in object code form under the terms
+of sections 4 and 5, provided that you also convey the
+machine-readable Corresponding Source under the terms of this License,
+in one of these ways:
+
+ a) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by the
+ Corresponding Source fixed on a durable physical medium
+ customarily used for software interchange.
+
+ b) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by a
+ written offer, valid for at least three years and valid for as
+ long as you offer spare parts or customer support for that product
+ model, to give anyone who possesses the object code either (1) a
+ copy of the Corresponding Source for all the software in the
+ product that is covered by this License, on a durable physical
+ medium customarily used for software interchange, for a price no
+ more than your reasonable cost of physically performing this
+ conveying of source, or (2) access to copy the
+ Corresponding Source from a network server at no charge.
+
+ c) Convey individual copies of the object code with a copy of the
+ written offer to provide the Corresponding Source. This
+ alternative is allowed only occasionally and noncommercially, and
+ only if you received the object code with such an offer, in accord
+ with subsection 6b.
+
+ d) Convey the object code by offering access from a designated
+ place (gratis or for a charge), and offer equivalent access to the
+ Corresponding Source in the same way through the same place at no
+ further charge. You need not require recipients to copy the
+ Corresponding Source along with the object code. If the place to
+ copy the object code is a network server, the Corresponding Source
+ may be on a different server (operated by you or a third party)
+ that supports equivalent copying facilities, provided you maintain
+ clear directions next to the object code saying where to find the
+ Corresponding Source. Regardless of what server hosts the
+ Corresponding Source, you remain obligated to ensure that it is
+ available for as long as needed to satisfy these requirements.
+
+ e) Convey the object code using peer-to-peer transmission, provided
+ you inform other peers where the object code and Corresponding
+ Source of the work are being offered to the general public at no
+ charge under subsection 6d.
+
+ A separable portion of the object code, whose source code is excluded
+from the Corresponding Source as a System Library, need not be
+included in conveying the object code work.
+
+ A "User Product" is either (1) a "consumer product", which means any
+tangible personal property which is normally used for personal, family,
+or household purposes, or (2) anything designed or sold for incorporation
+into a dwelling. In determining whether a product is a consumer product,
+doubtful cases shall be resolved in favor of coverage. For a particular
+product received by a particular user, "normally used" refers to a
+typical or common use of that class of product, regardless of the status
+of the particular user or of the way in which the particular user
+actually uses, or expects or is expected to use, the product. A product
+is a consumer product regardless of whether the product has substantial
+commercial, industrial or non-consumer uses, unless such uses represent
+the only significant mode of use of the product.
+
+ "Installation Information" for a User Product means any methods,
+procedures, authorization keys, or other information required to install
+and execute modified versions of a covered work in that User Product from
+a modified version of its Corresponding Source. The information must
+suffice to ensure that the continued functioning of the modified object
+code is in no case prevented or interfered with solely because
+modification has been made.
+
+ If you convey an object code work under this section in, or with, or
+specifically for use in, a User Product, and the conveying occurs as
+part of a transaction in which the right of possession and use of the
+User Product is transferred to the recipient in perpetuity or for a
+fixed term (regardless of how the transaction is characterized), the
+Corresponding Source conveyed under this section must be accompanied
+by the Installation Information. But this requirement does not apply
+if neither you nor any third party retains the ability to install
+modified object code on the User Product (for example, the work has
+been installed in ROM).
+
+ The requirement to provide Installation Information does not include a
+requirement to continue to provide support service, warranty, or updates
+for a work that has been modified or installed by the recipient, or for
+the User Product in which it has been modified or installed. Access to a
+network may be denied when the modification itself materially and
+adversely affects the operation of the network or violates the rules and
+protocols for communication across the network.
+
+ Corresponding Source conveyed, and Installation Information provided,
+in accord with this section must be in a format that is publicly
+documented (and with an implementation available to the public in
+source code form), and must require no special password or key for
+unpacking, reading or copying.
+
+ 7. Additional Terms.
+
+ "Additional permissions" are terms that supplement the terms of this
+License by making exceptions from one or more of its conditions.
+Additional permissions that are applicable to the entire Program shall
+be treated as though they were included in this License, to the extent
+that they are valid under applicable law. If additional permissions
+apply only to part of the Program, that part may be used separately
+under those permissions, but the entire Program remains governed by
+this License without regard to the additional permissions.
+
+ When you convey a copy of a covered work, you may at your option
+remove any additional permissions from that copy, or from any part of
+it. (Additional permissions may be written to require their own
+removal in certain cases when you modify the work.) You may place
+additional permissions on material, added by you to a covered work,
+for which you have or can give appropriate copyright permission.
+
+ Notwithstanding any other provision of this License, for material you
+add to a covered work, you may (if authorized by the copyright holders of
+that material) supplement the terms of this License with terms:
+
+ a) Disclaiming warranty or limiting liability differently from the
+ terms of sections 15 and 16 of this License; or
+
+ b) Requiring preservation of specified reasonable legal notices or
+ author attributions in that material or in the Appropriate Legal
+ Notices displayed by works containing it; or
+
+ c) Prohibiting misrepresentation of the origin of that material, or
+ requiring that modified versions of such material be marked in
+ reasonable ways as different from the original version; or
+
+ d) Limiting the use for publicity purposes of names of licensors or
+ authors of the material; or
+
+ e) Declining to grant rights under trademark law for use of some
+ trade names, trademarks, or service marks; or
+
+ f) Requiring indemnification of licensors and authors of that
+ material by anyone who conveys the material (or modified versions of
+ it) with contractual assumptions of liability to the recipient, for
+ any liability that these contractual assumptions directly impose on
+ those licensors and authors.
+
+ All other non-permissive additional terms are considered "further
+restrictions" within the meaning of section 10. If the Program as you
+received it, or any part of it, contains a notice stating that it is
+governed by this License along with a term that is a further
+restriction, you may remove that term. If a license document contains
+a further restriction but permits relicensing or conveying under this
+License, you may add to a covered work material governed by the terms
+of that license document, provided that the further restriction does
+not survive such relicensing or conveying.
+
+ If you add terms to a covered work in accord with this section, you
+must place, in the relevant source files, a statement of the
+additional terms that apply to those files, or a notice indicating
+where to find the applicable terms.
+
+ Additional terms, permissive or non-permissive, may be stated in the
+form of a separately written license, or stated as exceptions;
+the above requirements apply either way.
+
+ 8. Termination.
+
+ You may not propagate or modify a covered work except as expressly
+provided under this License. Any attempt otherwise to propagate or
+modify it is void, and will automatically terminate your rights under
+this License (including any patent licenses granted under the third
+paragraph of section 11).
+
+ However, if you cease all violation of this License, then your
+license from a particular copyright holder is reinstated (a)
+provisionally, unless and until the copyright holder explicitly and
+finally terminates your license, and (b) permanently, if the copyright
+holder fails to notify you of the violation by some reasonable means
+prior to 60 days after the cessation.
+
+ Moreover, your license from a particular copyright holder is
+reinstated permanently if the copyright holder notifies you of the
+violation by some reasonable means, this is the first time you have
+received notice of violation of this License (for any work) from that
+copyright holder, and you cure the violation prior to 30 days after
+your receipt of the notice.
+
+ Termination of your rights under this section does not terminate the
+licenses of parties who have received copies or rights from you under
+this License. If your rights have been terminated and not permanently
+reinstated, you do not qualify to receive new licenses for the same
+material under section 10.
+
+ 9. Acceptance Not Required for Having Copies.
+
+ You are not required to accept this License in order to receive or
+run a copy of the Program. Ancillary propagation of a covered work
+occurring solely as a consequence of using peer-to-peer transmission
+to receive a copy likewise does not require acceptance. However,
+nothing other than this License grants you permission to propagate or
+modify any covered work. These actions infringe copyright if you do
+not accept this License. Therefore, by modifying or propagating a
+covered work, you indicate your acceptance of this License to do so.
+
+ 10. Automatic Licensing of Downstream Recipients.
+
+ Each time you convey a covered work, the recipient automatically
+receives a license from the original licensors, to run, modify and
+propagate that work, subject to this License. You are not responsible
+for enforcing compliance by third parties with this License.
+
+ An "entity transaction" is a transaction transferring control of an
+organization, or substantially all assets of one, or subdividing an
+organization, or merging organizations. If propagation of a covered
+work results from an entity transaction, each party to that
+transaction who receives a copy of the work also receives whatever
+licenses to the work the party's predecessor in interest had or could
+give under the previous paragraph, plus a right to possession of the
+Corresponding Source of the work from the predecessor in interest, if
+the predecessor has it or can get it with reasonable efforts.
+
+ You may not impose any further restrictions on the exercise of the
+rights granted or affirmed under this License. For example, you may
+not impose a license fee, royalty, or other charge for exercise of
+rights granted under this License, and you may not initiate litigation
+(including a cross-claim or counterclaim in a lawsuit) alleging that
+any patent claim is infringed by making, using, selling, offering for
+sale, or importing the Program or any portion of it.
+
+ 11. Patents.
+
+ A "contributor" is a copyright holder who authorizes use under this
+License of the Program or a work on which the Program is based. The
+work thus licensed is called the contributor's "contributor version".
+
+ A contributor's "essential patent claims" are all patent claims
+owned or controlled by the contributor, whether already acquired or
+hereafter acquired, that would be infringed by some manner, permitted
+by this License, of making, using, or selling its contributor version,
+but do not include claims that would be infringed only as a
+consequence of further modification of the contributor version. For
+purposes of this definition, "control" includes the right to grant
+patent sublicenses in a manner consistent with the requirements of
+this License.
+
+ Each contributor grants you a non-exclusive, worldwide, royalty-free
+patent license under the contributor's essential patent claims, to
+make, use, sell, offer for sale, import and otherwise run, modify and
+propagate the contents of its contributor version.
+
+ In the following three paragraphs, a "patent license" is any express
+agreement or commitment, however denominated, not to enforce a patent
+(such as an express permission to practice a patent or covenant not to
+sue for patent infringement). To "grant" such a patent license to a
+party means to make such an agreement or commitment not to enforce a
+patent against the party.
+
+ If you convey a covered work, knowingly relying on a patent license,
+and the Corresponding Source of the work is not available for anyone
+to copy, free of charge and under the terms of this License, through a
+publicly available network server or other readily accessible means,
+then you must either (1) cause the Corresponding Source to be so
+available, or (2) arrange to deprive yourself of the benefit of the
+patent license for this particular work, or (3) arrange, in a manner
+consistent with the requirements of this License, to extend the patent
+license to downstream recipients. "Knowingly relying" means you have
+actual knowledge that, but for the patent license, your conveying the
+covered work in a country, or your recipient's use of the covered work
+in a country, would infringe one or more identifiable patents in that
+country that you have reason to believe are valid.
+
+ If, pursuant to or in connection with a single transaction or
+arrangement, you convey, or propagate by procuring conveyance of, a
+covered work, and grant a patent license to some of the parties
+receiving the covered work authorizing them to use, propagate, modify
+or convey a specific copy of the covered work, then the patent license
+you grant is automatically extended to all recipients of the covered
+work and works based on it.
+
+ A patent license is "discriminatory" if it does not include within
+the scope of its coverage, prohibits the exercise of, or is
+conditioned on the non-exercise of one or more of the rights that are
+specifically granted under this License. You may not convey a covered
+work if you are a party to an arrangement with a third party that is
+in the business of distributing software, under which you make payment
+to the third party based on the extent of your activity of conveying
+the work, and under which the third party grants, to any of the
+parties who would receive the covered work from you, a discriminatory
+patent license (a) in connection with copies of the covered work
+conveyed by you (or copies made from those copies), or (b) primarily
+for and in connection with specific products or compilations that
+contain the covered work, unless you entered into that arrangement,
+or that patent license was granted, prior to 28 March 2007.
+
+ Nothing in this License shall be construed as excluding or limiting
+any implied license or other defenses to infringement that may
+otherwise be available to you under applicable patent law.
+
+ 12. No Surrender of Others' Freedom.
+
+ If conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License. If you cannot convey a
+covered work so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you may
+not convey it at all. For example, if you agree to terms that obligate you
+to collect a royalty for further conveying from those to whom you convey
+the Program, the only way you could satisfy both those terms and this
+License would be to refrain entirely from conveying the Program.
+
+ 13. Use with the GNU Affero General Public License.
+
+ Notwithstanding any other provision of this License, you have
+permission to link or combine any covered work with a work licensed
+under version 3 of the GNU Affero General Public License into a single
+combined work, and to convey the resulting work. The terms of this
+License will continue to apply to the part which is the covered work,
+but the special requirements of the GNU Affero General Public License,
+section 13, concerning interaction through a network will apply to the
+combination as such.
+
+ 14. Revised Versions of this License.
+
+ The Free Software Foundation may publish revised and/or new versions of
+the GNU General Public License from time to time. Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+ Each version is given a distinguishing version number. If the
+Program specifies that a certain numbered version of the GNU General
+Public License "or any later version" applies to it, you have the
+option of following the terms and conditions either of that numbered
+version or of any later version published by the Free Software
+Foundation. If the Program does not specify a version number of the
+GNU General Public License, you may choose any version ever published
+by the Free Software Foundation.
+
+ If the Program specifies that a proxy can decide which future
+versions of the GNU General Public License can be used, that proxy's
+public statement of acceptance of a version permanently authorizes you
+to choose that version for the Program.
+
+ Later license versions may give you additional or different
+permissions. However, no additional obligations are imposed on any
+author or copyright holder as a result of your choosing to follow a
+later version.
+
+ 15. Disclaimer of Warranty.
+
+ THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
+APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
+IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+ 16. Limitation of Liability.
+
+ IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGES.
+
+ 17. Interpretation of Sections 15 and 16.
+
+ If the disclaimer of warranty and limitation of liability provided
+above cannot be given local legal effect according to their terms,
+reviewing courts shall apply local law that most closely approximates
+an absolute waiver of all civil liability in connection with the
+Program, unless a warranty or assumption of liability accompanies a
+copy of the Program in return for a fee.
+
+ END OF TERMS AND CONDITIONS
+
+ How to Apply These Terms to Your New Programs
+
+ If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+ To do so, attach the following notices to the program. It is safest
+to attach them to the start of each source file to most effectively
+state the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+
+ Copyright (C)
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see .
+
+Also add information on how to contact you by electronic and paper mail.
+
+ If the program does terminal interaction, make it output a short
+notice like this when it starts in an interactive mode:
+
+ Copyright (C)
+ This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+ This is free software, and you are welcome to redistribute it
+ under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License. Of course, your program's commands
+might be different; for a GUI interface, you would use an "about box".
+
+ You should also get your employer (if you work as a programmer) or school,
+if any, to sign a "copyright disclaimer" for the program, if necessary.
+For more information on this, and how to apply and follow the GNU GPL, see
+.
+
+ The GNU General Public License does not permit incorporating your program
+into proprietary programs. If your program is a subroutine library, you
+may consider it more useful to permit linking proprietary applications with
+the library. If this is what you want to do, use the GNU Lesser General
+Public License instead of this License. But first, please read
+.
+
diff --git a/src/external/CoMISo/CoMISo.cmake b/src/external/CoMISo/CoMISo.cmake
new file mode 100644
index 000000000..a981bf386
--- /dev/null
+++ b/src/external/CoMISo/CoMISo.cmake
@@ -0,0 +1,2 @@
+# add our macro directory to cmake search path
+set (CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}/cmake)
diff --git a/src/external/CoMISo/Config/CoMISoDefines.hh b/src/external/CoMISo/Config/CoMISoDefines.hh
new file mode 100644
index 000000000..492008b86
--- /dev/null
+++ b/src/external/CoMISo/Config/CoMISoDefines.hh
@@ -0,0 +1,48 @@
+/*===========================================================================*\
+ * *
+ * CoMISo *
+ * Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen *
+ * www.rwth-graphics.de *
+ * *
+ *---------------------------------------------------------------------------*
+ * This file is part of CoMISo. *
+ * *
+ * CoMISo is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * CoMISo is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with CoMISo. If not, see . *
+ * *
+\*===========================================================================*/
+
+#ifndef COMISODLLEXPORT
+ #ifdef WIN32
+ #ifdef COMISODLL
+ #ifdef USECOMISO
+ #define COMISODLLEXPORT __declspec(dllimport)
+ #define COMISODLLEXPORTONLY
+ #else
+ #define COMISODLLEXPORT __declspec(dllexport)
+ #define COMISODLLEXPORTONLY __declspec(dllexport)
+ #endif
+ #else
+ #define COMISODLLEXPORT
+ #define COMISODLLEXPORTONLY
+ #endif
+ #else
+ #define COMISODLLEXPORT
+ #define COMISODLLEXPORTONLY
+ #endif
+#endif
+
+#undef min
+#undef max
+
+
diff --git a/src/external/CoMISo/Config/config.hh.in b/src/external/CoMISo/Config/config.hh.in
new file mode 100644
index 000000000..fc1d67bf3
--- /dev/null
+++ b/src/external/CoMISo/Config/config.hh.in
@@ -0,0 +1,24 @@
+// Build time dependencies for CoMiso
+
+@COMISO_QT4_CONFIG_FILE_SETTINGS@
+@COMISO_BOOST_CONFIG_FILE_SETTINGS@
+@COMISO_BLAS_CONFIG_FILE_SETTINGS@
+@COMISO_GMM_CONFIG_FILE_SETTINGS@
+@COMISO_ADOLC_CONFIG_FILE_SETTINGS@
+@COMISO_SUITESPARSE_CONFIG_FILE_SETTINGS@
+@COMISO_SUITESPARSE_SPQR_CONFIG_FILE_SETTINGS@
+@COMISO_MPI_CONFIG_FILE_SETTINGS@
+@COMISO_HSL_CONFIG_FILE_SETTINGS@
+@COMISO_PETSC_CONFIG_FILE_SETTINGS@
+@COMISO_TAO_CONFIG_FILE_SETTINGS@
+@COMISO_IPOPT_CONFIG_FILE_SETTINGS@
+@COMISO_MUMPS_CONFIG_FILE_SETTINGS@
+@COMISO_METIS_CONFIG_FILE_SETTINGS@
+@COMISO_CGAL_CONFIG_FILE_SETTINGS@
+@COMISO_TAUCS_CONFIG_FILE_SETTINGS@
+@COMISO_GUROBI_CONFIG_FILE_SETTINGS@
+@COMISO_ARPACK_CONFIG_FILE_SETTINGS@
+@COMISO_CPLEX_CONFIG_FILE_SETTINGS@
+@COMISO_EIGEN3_CONFIG_FILE_SETTINGS@
+
+
diff --git a/src/external/CoMISo/EigenSolver/ArpackSolver.cc b/src/external/CoMISo/EigenSolver/ArpackSolver.cc
new file mode 100644
index 000000000..981ab2ffb
--- /dev/null
+++ b/src/external/CoMISo/EigenSolver/ArpackSolver.cc
@@ -0,0 +1,131 @@
+//=============================================================================
+//
+// CLASS ArpackSolver - IMPLEMENTATION
+//
+//=============================================================================
+
+#define COMISO_ARPACKSOLVER_C
+
+//== INCLUDES =================================================================
+
+#include "ArpackSolver.hh"
+
+//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
+#if (COMISO_ARPACK_AVAILABLE && COMISO_SUITESPARSE_AVAILABLE && COMISO_EIGEN3_AVAILABLE)
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== IMPLEMENTATION ==========================================================
+
+
+template
+void
+ArpackSolver::
+solve(const MatrixT& _A,
+ std::vector& _eigenvalues,
+ MatrixT2& _eigenvectors,
+ const int _n_eigvalues,
+ const char* _which_eigs )
+{
+ Matrix A(_A);
+// ARSymStdEig eig_prob(A.matrix().cols(), _n_eigvalues, &A, &Matrix::mult_Mv, (char*)_which_eigs,
+// 0, 0.0, 2000);
+
+ ARSymStdEig eig_prob(A.matrix().cols(), _n_eigvalues, &A, &Matrix::mult_Mv, (char*)_which_eigs,
+ 0, 0.0, 100000);
+
+ int n_converged = eig_prob.FindEigenvectors();
+
+ // store result
+ _eigenvalues.resize(n_converged);
+ _eigenvectors.resize(A.matrix().rows(),n_converged);
+
+ for( int i=0; i
+void
+ArpackSolver::
+solve_inverse(const MatrixT& _A,
+ std::vector& _eigenvalues,
+ MatrixT2& _eigenvectors,
+ const int _n_eigvalues,
+ const char* _which_eigs)
+{
+ Matrix A(_A,true);
+ ARSymStdEig eig_prob(A.matrix().cols(), _n_eigvalues, &A, &Matrix::mult_M_inv_v, (char*)_which_eigs,
+ 0, 0.0, 2000);
+
+ // ARSymStdEig(int np, int nevp, ARFOP* objOPp,
+ // void (ARFOP::* MultOPxp)(ARFLOAT[], ARFLOAT[]),
+ // char* whichp = "LM", int ncvp = 0, ARFLOAT tolp = 0.0,
+ // int maxitp = 0, ARFLOAT* residp = NULL, bool ishiftp = true);
+
+ int n_converged = eig_prob.FindEigenvectors();
+
+ // store result
+ _eigenvalues.resize(n_converged);
+ _eigenvectors.resize(A.matrix().rows(),n_converged);
+
+ for( int i=0; i
+void
+ArpackSolver::
+check_result(const MatrixT& _A, std::vector& _eigenvalues, MatrixT2& _eigenvectors)
+{
+ int n=_eigenvectors.rows();
+
+ if(n<20)
+ std::cerr << _A << std::endl;
+
+ for(unsigned int i=0; i<_eigenvalues.size(); ++i)
+ {
+ std::cerr << "eigenvalue " << i << ": " << _eigenvalues[i] << ", ";
+ if(n < 20)
+ {
+ std::cerr << "eigenvector: ";
+ for(int j=0; j v = _eigenvectors.block(0,i,n,1);
+ std::cerr << "residuum norm: " << (_A*v - _eigenvalues[i]*v).norm() << std::endl;
+ }
+}
+
+//-----------------------------------------------------------------------------
+
+
+
+//=============================================================================
+} // namespace COMISO
+//=============================================================================
+
+//=============================================================================
+#endif // COMISO_SUITESPARSE_AVAILABLE
+//=============================================================================
diff --git a/src/external/CoMISo/EigenSolver/ArpackSolver.hh b/src/external/CoMISo/EigenSolver/ArpackSolver.hh
new file mode 100644
index 000000000..067916a60
--- /dev/null
+++ b/src/external/CoMISo/EigenSolver/ArpackSolver.hh
@@ -0,0 +1,98 @@
+//=============================================================================
+//
+// CLASS ArpackSolver
+//
+//=============================================================================
+
+
+#ifndef COMISO_ARPACKSOLVER_HH
+#define COMISO_ARPACKSOLVER_HH
+
+//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
+#include
+#if (COMISO_ARPACK_AVAILABLE && COMISO_SUITESPARSE_AVAILABLE && COMISO_EIGEN3_AVAILABLE)
+
+//== INCLUDES =================================================================
+#include
+
+#include
+#include "EigenArpackMatrixT.hh"
+
+#include
+
+//== FORWARDDECLARATIONS ======================================================
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== CLASS DEFINITION =========================================================
+
+
+
+
+/** \class ArpackSolver ArpackSolver.hh
+
+ Brief Description.
+
+ A more elaborate description follows.
+*/
+
+
+class COMISODLLEXPORT ArpackSolver
+{
+public:
+
+ // sparse matrix type
+ typedef EigenArpackMatrixT > Matrix;
+
+
+ /// Constructor
+ ArpackSolver() {}
+
+ /// Destructor
+ ~ArpackSolver() {}
+
+ // solve eigenproblem
+ // number of desired eigenvalues -> _n_eigenvalues
+ // which eigenvalues -> one of {LA (largest algebraic), SA (smalles algebraic), LM (largest magnitude), SM(smallest magnitued), BE(both ends)}
+ template
+ void solve(const MatrixT& _A,
+ std::vector& _eigenvalues,
+ MatrixT2& _eigenvectors,
+ const int _n_eigvalues = 1,
+ const char* _which_eigs = "SM");
+
+ // solve eigenproblem
+ // number of desired eigenvalues -> _n_eigenvalues
+ // which eigenvalues -> one of {LA (largest algebraic), SA (smalles algebraic), LM (largest magnitude), SM(smallest magnitued), BE(both ends)}
+ template
+ void solve_inverse(const MatrixT& _A,
+ std::vector& _eigenvalues,
+ MatrixT2& _eigenvectors,
+ const int _n_eigvalues = 1,
+ const char* _which_eigs = "LM");
+
+
+ // check resulting eigenvalues/eigenvectors
+ template
+ void check_result(const MatrixT& _A, std::vector& _eigenvalues, MatrixT2& _eigenvectors);
+
+private:
+
+};
+
+
+//=============================================================================
+} // namespace ACG
+//=============================================================================
+#if defined(INCLUDE_TEMPLATES) && !defined(COMISO_ARPACKSOLVER_C)
+#define COMISO_ARPACKSOLVER_TEMPLATES
+#include "ArpackSolver.cc"
+#endif
+//=============================================================================
+#endif // COMISO_SUITESPARSE_AVAILABLE
+//=============================================================================
+#endif // ACG_ARPACKSOLVER_HH defined
+//=============================================================================
+
diff --git a/src/external/CoMISo/EigenSolver/EigenArpackMatrixT.cc b/src/external/CoMISo/EigenSolver/EigenArpackMatrixT.cc
new file mode 100644
index 000000000..540497817
--- /dev/null
+++ b/src/external/CoMISo/EigenSolver/EigenArpackMatrixT.cc
@@ -0,0 +1,35 @@
+//=============================================================================
+//
+// CLASS EigenArpackMatrixT - IMPLEMENTATION
+//
+//=============================================================================
+
+#define COMISO_EIGENARPACKMATRIXT_C
+
+//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
+#include
+#if COMISO_SUITESPARSE_AVAILABLE
+//=============================================================================
+
+
+//== INCLUDES =================================================================
+
+#include "EigenArpackMatrixT.hh"
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== IMPLEMENTATION ==========================================================
+
+
+
+//-----------------------------------------------------------------------------
+
+
+
+//=============================================================================
+} // namespace ACG
+//=============================================================================
+#endif // COMISO_SUITESPARSE_AVAILABLE
+//=============================================================================
diff --git a/src/external/CoMISo/EigenSolver/EigenArpackMatrixT.hh b/src/external/CoMISo/EigenSolver/EigenArpackMatrixT.hh
new file mode 100644
index 000000000..d7a112103
--- /dev/null
+++ b/src/external/CoMISo/EigenSolver/EigenArpackMatrixT.hh
@@ -0,0 +1,126 @@
+//=============================================================================
+//
+// CLASS EigenArpackMatrixT
+//
+//=============================================================================
+
+
+#ifndef COMISO_EIGENARPACKMATRIXT_HH
+#define COMISO_EIGENARPACKMATRIXT_HH
+
+//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
+#include
+#if (COMISO_SUITESPARSE_AVAILABLE && COMISO_EIGEN3_AVAILABLE)
+//=============================================================================
+
+
+//== INCLUDES =================================================================
+
+#include
+#include
+
+#if EIGEN_VERSION_AT_LEAST(3,1,0)
+ #include
+#else
+ #define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+ #include
+#endif
+#include
+
+//== FORWARDDECLARATIONS ======================================================
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== CLASS DEFINITION =========================================================
+
+
+/** \class EigenArpackMatrixT EigenArpackMatrixT.hh
+
+ Brief Description.
+
+ A more elaborate description follows.
+*/
+
+template
+class EigenArpackMatrixT
+{
+public:
+
+ typedef MatrixT Matrix;
+ typedef RealT Real;
+
+ /// Constructor
+ template
+ EigenArpackMatrixT(const MatrixT2& _m, bool _use_inverse = false)
+ {
+ mat_ = _m;
+
+ if(_use_inverse)
+ {
+ sllt_.compute(mat_);
+
+#if EIGEN_VERSION_AT_LEAST(3,1,0)
+ if ( !sllt_.info() != Eigen::Success )
+#else
+ if ( !sllt_.succeeded() )
+#endif
+ std::cout << "[ERROR] EigenArpackMatrix(): Could not compute llt factorization." << std::endl;
+ }
+ }
+
+ /// Destructor
+ ~EigenArpackMatrixT() {}
+
+ // get reference on matrix
+ Matrix& matrix() { return mat_; }
+
+ // matrix-vector multiplication _w = mat_*_v
+ void mult_Mv(Real* _v, Real* _w)
+ {
+ Eigen::Map > v(_v,mat_.rows()); // uses v as a ArrayXf object
+ Eigen::Map > w(_w,mat_.cols()); // uses w as a ArrayXf object
+
+ w = mat_*v;
+ }
+
+ // matrix-vector multiplication _w = mat_*_v
+ void mult_M_inv_v(Real* _v, Real* _w)
+ {
+ Eigen::Map > v(_v,mat_.rows()); // uses v as a ArrayXf object
+ Eigen::Map > w(_w,mat_.cols()); // uses w as a ArrayXf object
+
+ w = sllt_.solve(v);
+
+// std::cerr << "input:" << std::endl;
+// std::cerr << v << std::endl;
+// std::cerr << "output:" << std::endl;
+// std::cerr << w << std::endl;
+ }
+
+private:
+
+ Matrix mat_;
+
+#if EIGEN_VERSION_AT_LEAST(3,1,0)
+ Eigen::CholmodSupernodalLLT > sllt_;
+#else
+ Eigen::SparseLLT, Eigen::Cholmod> sllt_;
+#endif
+};
+
+
+//=============================================================================
+} // namespace COMISO
+//=============================================================================
+#if defined(INCLUDE_TEMPLATES) && !defined(COMISO_EIGENARPACKMATRIXT_C)
+#define COMISO_EIGENARPACKMATRIXT_TEMPLATES
+#include "EigenArpackMatrixT.cc"
+#endif
+//=============================================================================
+#endif // COMISO_SUITESPARSE_AVAILABLE
+//=============================================================================
+#endif // COMISO_EIGENARPACKMATRIXT_HH defined
+//=============================================================================
+
diff --git a/src/external/CoMISo/Examples/factored_solver/CMakeLists.txt b/src/external/CoMISo/Examples/factored_solver/CMakeLists.txt
new file mode 100644
index 000000000..da380df66
--- /dev/null
+++ b/src/external/CoMISo/Examples/factored_solver/CMakeLists.txt
@@ -0,0 +1,15 @@
+include (CoMISoExample)
+
+if (WIN32)
+ acg_add_executable (factored_solver WIN32 ${sources} ${headers} )
+else ()
+ acg_add_executable (factored_solver ${sources} ${headers} )
+endif ()
+
+# enable rpath linking
+set_target_properties(factored_solver PROPERTIES INSTALL_RPATH_USE_LINK_PATH 1)
+
+target_link_libraries (factored_solver
+ CoMISo
+ ${COMISO_LINK_LIBRARIES}
+)
diff --git a/src/external/CoMISo/Examples/factored_solver/main.cc b/src/external/CoMISo/Examples/factored_solver/main.cc
new file mode 100644
index 000000000..9248aac6f
--- /dev/null
+++ b/src/external/CoMISo/Examples/factored_solver/main.cc
@@ -0,0 +1,193 @@
+/*===========================================================================*\
+ * *
+ * CoMISo *
+ * Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen *
+ * www.rwth-graphics.de *
+ * *
+ *---------------------------------------------------------------------------*
+ * This file is part of CoMISo. *
+ * *
+ * CoMISo is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * CoMISo is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with CoMISo. If not, see . *
+ * *
+\*===========================================================================*/
+
+
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+/// function to setup a random sparse row matrix of dimension _m x _n
+/// for the simplicity of this example only integer valued entries are used
+template
+void random_sparse_row_matrix( MatrixT& _B, int _m, int _n, double _density = 0.7)
+{
+ gmm::resize(_B, _m, _n);
+
+ for( int i=0; i<_m; ++i)
+ for( int j=0; j<_n; ++j)
+ if( (rand()-1.0*_density*RAND_MAX)/RAND_MAX> 0) // for sparseness
+ _B(i,j) = round(((rand()-0.4*RAND_MAX)/RAND_MAX)*10.0);
+}
+
+/// function to setup a random sparse constraint row matrix of dimension _c x _n
+/// for the simplicity of the example only -1, 0, 1 constraints are used
+template
+void simple_constraint_row_matrix( MatrixT& _C, int _c, int _n, double _distribution = 0.2)
+{
+ gmm::resize( _C, _c, _n);
+ for( int i=0; i<_c; ++i)
+ for( int j=0; j<_n; ++j)
+ {
+ double randnum = (double(rand())/double(RAND_MAX));
+ if ( randnum < _distribution)
+ _C( i,j) = -1;
+ else if( randnum > (1.0-_distribution))
+ _C( i,j) = 1;
+ else
+ _C( i,j) = 0;
+ }
+}
+
+/// function to print the equations corresponding to the matrices of an equation system
+template
+void print_equations( const MatrixT& _B)
+{
+ int m = gmm::mat_nrows( _B);
+ int n = gmm::mat_ncols( _B);
+ for( int i = 0; i < m; ++i)
+ {
+ for( int j = 0; j < n-1; ++j)
+ {
+ if( _B(i,j) != 0.0)
+ std::cout << _B(i,j) << "*x" << j;
+ else
+ std::cout << " 0 ";
+ if( j < n-2 ) std::cout << " + ";
+ }
+ std::cout << " = " << _B(i, n-1) << std::endl;
+ }
+}
+
+
+// Example main
+int main(void)
+{
+ std::cout << "---------- 1) setup an (m x n) sparse row matrix B (i.e. the B in the system ((Bx)^T)Bx)" << std::endl;
+ int m = 9;
+ int n = 5+1;
+ gmm::row_matrix< gmm::wsvector< double > > B;
+ random_sparse_row_matrix( B, m, n, 0.85);
+ std::cout << B << std::endl << std::endl;
+ //gmm::inspect_matrix( B );
+ std::cout << std::endl;
+
+ std::cout << "---------- 2) define a set of linear constraints as an (c x n) row matrix C" << std::endl;
+ int c = 2;
+ gmm::row_matrix< gmm::wsvector< double > > C;
+ simple_constraint_row_matrix( C, c, n);
+ std::cout << C << std::endl;
+ std::cout << "corresponding to the following linear equations : " << std::endl;
+ print_equations( C );
+ std::cout << std::endl;
+
+ std::cout << "---------- 3) we now explicitly carry out the steps performed internally by the constrained solver and compare the two results at the end..." << std::endl;
+ // copy the matrices
+ gmm::row_matrix< gmm::wsvector< double > > Bcpy( B );
+ gmm::row_matrix< gmm::wsvector< double > > Ccpy( C );
+
+ // create a constrained solver
+ COMISO::ConstrainedSolver cs;
+ // vector of indices to round (this is the mixed-integer part)
+ std::vector< int > ids_to_round;
+ // lets say we want to round the third variable
+ ids_to_round.push_back(2);
+
+ // vector of independent variables to be eliminated (computed by the make_constraints_independent function)
+ std::vector< int > ids_to_elim;
+
+ std::cout << "---------- ---------- 3.1) make the constraints independent (gauss elimination on C)" << std::endl;
+ print_equations( Ccpy );
+ cs.make_constraints_independent( Ccpy, ids_to_round, ids_to_elim);
+ std::cout << " constraints after gauss elimination..." << std::endl;
+ std::cout << Ccpy << std::endl;
+ std::cout << " the variables to be eliminated are: " << std::endl;
+ std::cout << ids_to_elim << std::endl << std::endl;
+
+ std::cout << "---------- ---------- 3.2) eliminate constraints from system matrix B" << std::endl;
+ // this is the column matrix later used by the solver, it is setup by the eliminate_constraints function
+ gmm::col_matrix< gmm::wsvector< double > > Bcol;
+
+ // this re-indexing is also used by the solver, to know which variables are still there (!=-1) and which have been eliminated (=-1) it is setup by eliminate_constraints
+ std::vector< int > new_idx;
+
+ cs.eliminate_constraints( Ccpy, Bcpy, ids_to_round, ids_to_elim, new_idx, Bcol);
+ std::cout << " B matrix after elimination of constraints..." << std::endl;
+ std::cout << Bcol << std::endl;
+
+ std::cout << "---------- ---------- 3.3) setup the linear system Ax=b, where by forming B^TB and extracting the right hand side" << std::endl;
+
+ // this is the solution vector x
+ std::vector< double > x;
+
+ int new_n = gmm::mat_ncols( Bcol);
+ // set up B transposed
+ gmm::col_matrix< gmm::wsvector< double > > Bt( new_n, m);
+ gmm::copy( gmm::transposed( Bcol), Bt);
+
+ // setup BtB
+ gmm::col_matrix< gmm::wsvector< double > > BtB( new_n, new_n);
+ gmm::mult( Bt, Bcol, BtB);
+
+ // extract rhs
+ std::vector< double > rhs( new_n);
+ gmm::copy( gmm::scaled(gmm::mat_const_col( BtB, new_n - 1),-1.0), rhs);
+ rhs.resize( new_n - 1);
+
+ // resize BtB to only contain the actual system matrix (and not the rhs)
+ gmm::resize( BtB, new_n - 1, new_n - 1);
+ x.resize( new_n - 1);
+
+ // BtB -> CSC
+ gmm::csc_matrix BtBCSC;
+ BtBCSC.init_with_good_format( BtB);
+
+ std::cout << " the linear system now looks like..." << std::endl;
+ std::cout << " Matrix A\n " << BtBCSC << std::endl;
+ std::cout << " Right hand side b\n" << rhs << std::endl << std::endl;
+
+ std::cout << "---------- ---------- 3.4) solve the system using the mixed-integer solver..." << std::endl;
+ // create solver
+ COMISO::MISolver miso;
+ // miso solve
+ miso.solve( BtBCSC, x, rhs, ids_to_round);
+ std::cout << " solution vector x is\n" << x << std::endl << std::endl;
+
+ std::cout << "---------- ---------- 3.5) now the solution must be re-indexed to the expected/original form/size...." << std::endl;
+ cs.restore_eliminated_vars( Ccpy, x, ids_to_elim, new_idx);
+ std::cout << " fullsize solution vector x is\n" << x << std::endl << std::endl;
+
+ std::cout << "---------- ---------- 4) the same result is obtained by one call to the constrained solver, which takes care of re-indexing etc. internally..." << std::endl;
+ // ids_to_round is altered by previous steps...
+ ids_to_round.clear();
+ ids_to_round.push_back(2);
+ cs.solve( C, B, x, ids_to_round, 0.0, false, true);
+ std::cout << " fullsize solution vector x is\n" << x << std::endl << std::endl;
+
+ return -1;
+}
+
diff --git a/src/external/CoMISo/Examples/quadratic_solver/CMakeLists.txt b/src/external/CoMISo/Examples/quadratic_solver/CMakeLists.txt
new file mode 100644
index 000000000..a76e75913
--- /dev/null
+++ b/src/external/CoMISo/Examples/quadratic_solver/CMakeLists.txt
@@ -0,0 +1,15 @@
+include (CoMISoExample)
+
+if (WIN32)
+ acg_add_executable (quadratic_solver WIN32 ${sources} ${headers} )
+else ()
+ acg_add_executable (quadratic_solver ${sources} ${headers} )
+endif ()
+
+# enable rpath linking
+set_target_properties(quadratic_solver PROPERTIES INSTALL_RPATH_USE_LINK_PATH 1)
+
+target_link_libraries (quadratic_solver
+ CoMISo
+ ${COMISO_LINK_LIBRARIES}
+)
diff --git a/src/external/CoMISo/Examples/quadratic_solver/main.cc b/src/external/CoMISo/Examples/quadratic_solver/main.cc
new file mode 100644
index 000000000..74675be6d
--- /dev/null
+++ b/src/external/CoMISo/Examples/quadratic_solver/main.cc
@@ -0,0 +1,199 @@
+/*===========================================================================*\
+ * *
+ * CoMISo *
+ * Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen *
+ * www.rwth-graphics.de *
+ * *
+ *---------------------------------------------------------------------------*
+ * This file is part of CoMISo. *
+ * *
+ * CoMISo is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * CoMISo is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with CoMISo. If not, see . *
+ * *
+\*===========================================================================*/
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+/// function to setup a random sparse row matrix of dimension _m x _n
+/// for the simplicity of this example only integer valued entries are used
+template
+void random_sparse_row_matrix( MatrixT& _B, int _m, int _n, double _density = 0.7)
+{
+ gmm::resize(_B, _m, _n);
+
+ for( int i=0; i<_m; ++i)
+ for( int j=0; j<_n; ++j)
+ if( (rand()-1.0*_density*RAND_MAX)/RAND_MAX> 0) // for sparseness
+ _B(i,j) = round(((rand()-0.4*RAND_MAX)/RAND_MAX)*10.0);
+}
+
+/// function to extract the actual system Ax=b of linear equation from a B^tB matrix
+template
+void extract_Axb( const RMatrixT& _B, CMatrixT& _A, std::vector< double >& _b)
+{
+ int dimm = gmm::mat_nrows(_B);
+ int dimn = gmm::mat_ncols(_B);
+ gmm::col_matrix< gmm::wsvector< double > > Btcol;
+ gmm::col_matrix< gmm::wsvector< double > > Bcol;
+ gmm::resize( Btcol, dimn, dimm);
+ gmm::resize( Bcol, dimm, dimn);
+ gmm::resize( _A, dimn, dimn);
+ gmm::copy( _B, Bcol);
+ gmm::copy( gmm::transposed( Bcol), Btcol);
+ gmm::mult( Btcol, Bcol, _A);
+ _b.resize( dimn);
+ gmm::copy( _A.col(dimn-1), _b);
+ _b.resize( dimn-1);
+ gmm::resize( _A, dimn-1, dimn-1);
+ gmm::scale(_b, -1.0);
+}
+
+/// function to setup a random sparse constraint row matrix of dimension _c x _n
+/// for the simplicity of the example only -1, 0, 1 constraints are used
+template
+void simple_constraint_row_matrix( MatrixT& _C, int _c, int _n, double _distribution = 0.2)
+{
+ gmm::resize( _C, _c, _n);
+ for( int i=0; i<_c; ++i)
+ for( int j=0; j<_n; ++j)
+ {
+ double randnum = (double(rand())/double(RAND_MAX));
+ if ( randnum < _distribution)
+ _C( i,j) = -1;
+ else if( randnum > (1.0-_distribution))
+ _C( i,j) = 1;
+ else
+ _C( i,j) = 0;
+ }
+}
+
+/// function to print the equations corresponding to the matrices of an equation system
+template
+void print_equations( const MatrixT& _B)
+{
+ int m = gmm::mat_nrows( _B);
+ int n = gmm::mat_ncols( _B);
+ for( int i = 0; i < m; ++i)
+ {
+ for( int j = 0; j < n-1; ++j)
+ {
+ if( _B(i,j) != 0.0)
+ std::cout << _B(i,j) << "*x" << j;
+ else
+ std::cout << " 0 ";
+ if( j < n-2 ) std::cout << " + ";
+ }
+ std::cout << " = " << _B(i, n-1) << std::endl;
+ }
+}
+
+
+// Example main
+int main(void)
+{
+ std::cout << "---------- 1) setup an (m x n) sparse row matrix B (i.e. the B in the system ((Bx)^T)Bx)" << std::endl;
+ int m = 9;
+ int n = 5+1;
+ gmm::row_matrix< gmm::wsvector< double > > B;
+ random_sparse_row_matrix( B, m, n, 0.85);
+ std::cout << B << std::endl << std::endl;
+
+ std::cout << "---------- 2) extract the Ax=b equation system, A (n-1 x n-1)" << std::endl;
+ gmm::col_matrix< gmm::wsvector< double > > A;
+ std::vector< double > b;
+ extract_Axb( B, A, b);
+ std::cout << " A " << std::endl;
+ std::cout << A << " " << b << std::endl;
+
+ //gmm::inspect_matrix( B );
+ std::cout << std::endl;
+
+ std::cout << "---------- 3) define a set of linear constraints as an (c x n) row matrix C" << std::endl;
+ int c = 2;
+ gmm::row_matrix< gmm::wsvector< double > > C;
+ simple_constraint_row_matrix( C, c, n);
+ std::cout << C << std::endl;
+ std::cout << "corresponding to the following linear equations : " << std::endl;
+ print_equations( C );
+ std::cout << std::endl;
+
+ std::cout << "---------- 4) we now explicitly carry out the steps performed internally by the constrained solver and compare the two results at the end..." << std::endl;
+ // copy the matrices
+ gmm::col_matrix< gmm::wsvector< double > > Acpy( A );
+
+ // create a constrained solver
+ COMISO::ConstrainedSolver cs;
+ // vector of indices to round (this is the mixed-integer part)
+ std::vector< int > ids_to_round;
+ // lets say we want to round the third variable
+ ids_to_round.push_back(2);
+
+ // vector of independent variables to be eliminated (computed by the make_constraints_independent function)
+ std::vector< int > ids_to_elim;
+
+ std::cout << "---------- ---------- 4.1) make the constraints independent (gauss elimination on C)" << std::endl;
+ print_equations( C );
+ cs.make_constraints_independent( C, ids_to_round, ids_to_elim);
+ std::cout << " constraints after gauss elimination..." << std::endl;
+ std::cout << C << std::endl;
+ std::cout << " the variables to be eliminated are: " << std::endl;
+ std::cout << ids_to_elim << std::endl << std::endl;
+ gmm::row_matrix< gmm::wsvector< double > > Ccpy( C );
+
+
+ std::cout << "---------- ---------- 4.2) eliminate constraints from system matrix A" << std::endl;
+
+ // CSC matrix later initialized and used by solver
+ gmm::csc_matrix< double > Acsc;
+
+ // this re-indexing is also used by the solver, to know which variables are still there (!=-1) and which have been eliminated (=-1) it is setup by eliminate_constraints
+ std::vector< int > new_idx;
+
+ std::vector< double > x(b.size());
+ std::vector< double > b_cpy(b);
+
+ cs.eliminate_constraints( Ccpy, Acpy, x, b, ids_to_round, ids_to_elim, new_idx, Acsc);
+ std::cout << " A matrix after elimination of constraints..." << std::endl;
+ std::cout << Acsc << std::endl;
+
+
+ std::cout << "---------- ---------- 4.3) solve the system using the mixed-integer solver..." << std::endl;
+ // create solver
+ COMISO::MISolver miso;
+ // miso solve
+ miso.solve( Acsc, x, b, ids_to_round);
+ std::cout << " solution vector x is\n" << x << std::endl << std::endl;
+
+
+ std::cout << "---------- ---------- 4.4) now the solution must be re-indexed to the expected/original form/size...." << std::endl;
+ cs.restore_eliminated_vars( Ccpy, x, ids_to_elim, new_idx);
+ std::cout << " fullsize solution vector x is\n" << x << std::endl << std::endl;
+
+
+ std::cout << "---------- ---------- 5) the same result is obtained by one call to the constrained solver, which takes care of re-indexing etc. internally..." << std::endl;
+ // ids_to_round is altered by previous steps...
+ ids_to_round.clear();
+ ids_to_round.push_back(2);
+ x.resize(gmm::mat_nrows(A));
+ b.resize(gmm::mat_nrows(A));
+ cs.solve( C, A, x, b_cpy, ids_to_round, 0.0, false, true);
+ std::cout << " fullsize solution vector x is\n" << x << std::endl << std::endl;
+
+ return -1;
+}
+
diff --git a/src/external/CoMISo/Examples/small_cplex_soc/CMakeLists.txt b/src/external/CoMISo/Examples/small_cplex_soc/CMakeLists.txt
new file mode 100644
index 000000000..b76596e7e
--- /dev/null
+++ b/src/external/CoMISo/Examples/small_cplex_soc/CMakeLists.txt
@@ -0,0 +1,15 @@
+include (CoMISoExample)
+
+if (WIN32)
+ acg_add_executable (small_cplex_soc WIN32 ${sources} ${headers} )
+else ()
+ acg_add_executable (small_cplex_soc ${sources} ${headers} )
+endif ()
+
+# enable rpath linking
+set_target_properties(small_cplex_soc PROPERTIES INSTALL_RPATH_USE_LINK_PATH 1)
+
+target_link_libraries (small_cplex_soc
+ CoMISo
+ ${COMISO_LINK_LIBRARIES}
+)
diff --git a/src/external/CoMISo/Examples/small_cplex_soc/main.cc b/src/external/CoMISo/Examples/small_cplex_soc/main.cc
new file mode 100644
index 000000000..53625f53e
--- /dev/null
+++ b/src/external/CoMISo/Examples/small_cplex_soc/main.cc
@@ -0,0 +1,124 @@
+/*===========================================================================*\
+ * *
+ * CoMISo *
+ * Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen *
+ * www.rwth-graphics.de *
+ * *
+ *---------------------------------------------------------------------------*
+ * This file is part of CoMISo. *
+ * *
+ * CoMISo is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * CoMISo is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with CoMISo. If not, see . *
+ * *
+\*===========================================================================*/
+
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+// solve least squares problem for x=1, y=2 and x-2y+z = 1
+// with hard constraints x =-3, z>=3, z^2 >= x^2+y^2
+
+
+// Example main
+int main(void)
+{
+ std::cout << "---------- 1) Problem description..." << std::endl;
+ std::cout << "Least squares terms: x=1, y=2 and x-2y+z = 1" << std::endl;
+ std::cout << "Constraints : x =-3, z>=3, z^2 >= x^2+y^2" << std::endl;
+
+ std::cout << "---------- 1) Get an instance of a LeastSquaresProblem..." << std::endl;
+ // number of unknowns
+ const int n = 3;
+ COMISO::LeastSquaresProblem lsqp(n);
+
+ // term0
+ COMISO::LinearConstraint::SVectorNC coeffs0(n);
+ coeffs0.coeffRef(0) = 1.0;
+ COMISO::LinearConstraint term0(coeffs0,-1.0,COMISO::NConstraintInterface::NC_EQUAL);
+ lsqp.add_term(&term0);
+
+ // term1
+ COMISO::LinearConstraint::SVectorNC coeffs1(n);
+ coeffs1.coeffRef(1) = 1.0;
+ COMISO::LinearConstraint term1(coeffs1,-2.0,COMISO::NConstraintInterface::NC_EQUAL);
+ lsqp.add_term(&term1);
+
+ // term2
+ COMISO::LinearConstraint::SVectorNC coeffs2(n);
+ coeffs2.coeffRef(0) = 1.0;
+ coeffs2.coeffRef(1) = -2.0;
+ coeffs2.coeffRef(2) = 1.0;
+ COMISO::LinearConstraint term2(coeffs2,-1.0,COMISO::NConstraintInterface::NC_EQUAL);
+ lsqp.add_term(&term2);
+
+ std::cout << "---------- 2) set up constraints" << std::endl;
+
+ // set x = -3.0
+ COMISO::LinearConstraint lc;
+ lc.coeffs().coeffRef(0) = 1.0;
+ lc.b() = 3.0;
+
+ // set z>=3 (cone constraint requires that z>=0 !!!)
+ COMISO::BoundConstraint bc(2,3,3,COMISO::NConstraintInterface::NC_GREATER_EQUAL);
+
+ // set z^2 >= x^2+y^2
+ COMISO::ConeConstraint cc;
+ cc.resize(3);
+ cc.i() = 2;
+ cc.c() = 4.0;
+ cc.Q()(0,0) = 2.0;
+ cc.Q()(1,1) = 4.0;
+ cc.Q()(0,1) = 1.0;
+ cc.Q()(1,0) = 1.0;
+
+ // fill constraint vector
+ std::vector constraints;
+ constraints.push_back(&lc);
+ constraints.push_back(&bc);
+ constraints.push_back(&cc);
+
+// check if CPLEX solver available in current configuration
+#if( COMISO_CPLEX_AVAILABLE)
+ std::cout << "---------- 3) Solve with CPLEX solver... " << std::endl;
+
+ COMISO::CPLEXSolver cplx;
+ cplx.solve(&lsqp, constraints);
+#endif
+
+ std::cout << "---------- 4) Print solution CPLEX..." << std::endl;
+ for( int i=0; iSetStringValue("derivative_test", "second-order");
+ ipopt.solve(&lsqp, constraints);
+#endif
+
+ std::cout << "---------- 6) Print solution..." << std::endl;
+ for( int i=0; i 0)
+ {
+ A.coeffRef(i,i-1) = -1.0;
+ ++count;
+ }
+ if(i evals;
+ Matrix evects;
+ arsolv.solve(A, evals, evects, m);
+
+ std::cout << "---------- 3) printing results..." << std::endl;
+ std::cerr << "********* eigenvalues: ";
+ for(unsigned int i=0; i. *
+ * *
+\*===========================================================================*/
+
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+
+/// function to initialize a simple row matrix of equations
+template
+void init_fac( MatrixT& _B )
+{
+ _B(0,0) = 0 ; _B(0,1) = 4 ; _B(0,2) = -2 ; _B(0,3) = 0 ; _B(0,4) = -1;
+ _B(1,0) = 0 ; _B(1,1) = 0 ; _B(1,2) = 0 ; _B(1,3) = 0 ; _B(1,4) = 0 ;
+ _B(2,0) = 5 ; _B(2,1) = 0 ; _B(2,2) = -3 ; _B(2,3) = 0 ; _B(2,4) = 0 ;
+ _B(3,0) = 0 ; _B(3,1) = 0 ; _B(3,2) = -2 ; _B(3,3) = 0 ; _B(3,4) = 0 ;
+ _B(4,0) = 0 ; _B(4,1) = -2; _B(4,2) = 0 ; _B(4,3) = 2 ; _B(4,4) = 0 ;
+}
+
+/// function to print the equations corresponding to the matrices of an equation system
+template
+void print_equations( const MatrixT& _B)
+{
+ int m = gmm::mat_nrows( _B);
+ int n = gmm::mat_ncols( _B);
+ for( int i = 0; i < m; ++i)
+ {
+ for( int j = 0; j < n-1; ++j)
+ {
+ if( _B(i,j) != 0.0)
+ std::cout << _B(i,j) << "*x" << j;
+ else
+ std::cout << " 0 ";
+ if( j < n-2 ) std::cout << " + ";
+ }
+ std::cout << " = " << _B(i, n-1) << std::endl;
+ }
+}
+
+
+// Example main
+int main(void)
+{
+ std::cout << "---------- 1) Setup a number of equations (i.e. the B matrix of a factored system of linear equations B^tB)..." << std::endl;
+ int n = 4;
+ gmm::row_matrix< gmm::wsvector< double > > B(n+1,n+1);
+ std::vector< double > x(n);
+
+ init_fac( B );
+
+ // create an empty constraint matrix (will be used later)
+ gmm::row_matrix< gmm::wsvector< double > > constraints(0,n+1); //n+1 because of right hand side
+ // create an empty vector of variable indices to be rounded (will be used later)
+ std::vector< int > ids_to_round;
+
+ std::cout << B << std::endl << std::endl;
+
+
+ std::cout << "---------- 2) The original solution to this system is..." << std::endl;
+
+ COMISO::ConstrainedSolver cs;
+ // void solve( RMatrixT& _constraints, RMatrixT& _B, VectorT& _x, VectorIT& _idx_to_round, double _reg_factor = 0.0, bool _show_miso_settings = true, bool _show_timings = true );
+ //_show_miso_settings requires a QT context and hence must be false in this example
+ cs.solve( constraints, B, x, ids_to_round, 0.0, false, true);
+ // copy this solution for later
+ std::vector< double > org_x( x);
+ std::cout << x << std::endl;
+
+
+ std::cout << "---------- 3) Rounding: forcing the second variable to lie on an integer, changes the solution to..." << std::endl;
+ // reset system
+ init_fac( B );
+ ids_to_round.push_back(1);
+ cs.solve( constraints, B, x, ids_to_round, 0.0, false, true);
+ std::cout << x << std::endl;
+
+
+ std::cout << "---------- 4) Constraining: forcing the first variable to equal the second changes the solution to..." << std::endl;
+ // reset system
+ init_fac( B );
+ ids_to_round.clear();
+ ids_to_round.push_back(1);
+ // setup constraint x0*1+x1*0+x2*(-1)+x3*0=0
+ gmm::resize( constraints, 1, n+1);
+ constraints( 0, 0 ) = 1.0;
+ constraints( 0, 1 ) = -1.0;
+ std::cout << " the constraint equation looks like this:" << std::endl;
+ print_equations( constraints);
+ cs.solve( constraints, B, x, ids_to_round, 0.0, false, true);
+ std::cout << x << std::endl;
+
+ return -1;
+}
+
diff --git a/src/external/CoMISo/Examples/small_miqp/CMakeLists.txt b/src/external/CoMISo/Examples/small_miqp/CMakeLists.txt
new file mode 100644
index 000000000..f4040ef8e
--- /dev/null
+++ b/src/external/CoMISo/Examples/small_miqp/CMakeLists.txt
@@ -0,0 +1,15 @@
+include (CoMISoExample)
+
+if (WIN32)
+ acg_add_executable (small_miqp WIN32 ${sources} ${headers} )
+else ()
+ acg_add_executable (small_miqp ${sources} ${headers} )
+endif ()
+
+# enable rpath linking
+set_target_properties(small_miqp PROPERTIES INSTALL_RPATH_USE_LINK_PATH 1)
+
+target_link_libraries (small_miqp
+ CoMISo
+ ${COMISO_LINK_LIBRARIES}
+)
diff --git a/src/external/CoMISo/Examples/small_miqp/main.cc b/src/external/CoMISo/Examples/small_miqp/main.cc
new file mode 100644
index 000000000..54a07e7cf
--- /dev/null
+++ b/src/external/CoMISo/Examples/small_miqp/main.cc
@@ -0,0 +1,156 @@
+/*===========================================================================*\
+ * *
+ * CoMISo *
+ * Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen *
+ * www.rwth-graphics.de *
+ * *
+ *---------------------------------------------------------------------------*
+ * This file is part of CoMISo. *
+ * *
+ * CoMISo is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * CoMISo is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with CoMISo. If not, see . *
+ * *
+\*===========================================================================*/
+
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+
+// generate an instance of a nonlinear problem by deriving from base class NProblemInterface
+// implement all virtual functions in order to solve this problem by any of the solvers located
+// in CoMISo/NSolver
+
+class SmallNProblem : public COMISO::NProblemInterface
+{
+public:
+
+ // Sparse Matrix Type
+ // typedef Eigen::DynamicSparseMatrix SMatrixNP;
+
+
+ // specify a function which has several local minima
+ // f(x,y)=(x-2y+1)^2 + (x-5)^2
+
+ // number of unknown variables, here x and y = 2
+ virtual int n_unknowns ( )
+ {
+ return 2;
+ }
+
+ // initial value where the optimization should start from
+ virtual void initial_x ( double* _x )
+ {
+ _x[0] = 0.0;
+ _x[1] = 0.0;
+ }
+
+ // function evaluation at location _x
+ virtual double eval_f ( const double* _x )
+ {
+ double term = _x[0] - 2.0*_x[1] + 1.0;
+ double term2 = _x[0] - 5.0;
+
+ return term*term + term2*term2;
+ }
+
+ // gradient evaluation at location _x
+ virtual void eval_gradient( const double* _x, double* _g)
+ {
+ double term = _x[0] - 2.0*_x[1] + 1.0;
+ double term2 = _x[0] - 5.0;
+
+ _g[0] = 2.0*term + 2.0*term2;
+ _g[1] = -4.0*term;
+ }
+
+ // hessian matrix evaluation at location _x
+ virtual void eval_hessian ( const double* _x, SMatrixNP& _H)
+ {
+ _H.resize(n_unknowns(), n_unknowns());
+ _H.setZero();
+
+ _H.coeffRef(0,0) = 4.0;
+ _H.coeffRef(1,0) = -4.0;
+ _H.coeffRef(0,1) = -4.0;
+ _H.coeffRef(1,1) = 8.0;
+ }
+
+ // print result
+ virtual void store_result ( const double* _x )
+ {
+ std::cerr << "Energy: " << eval_f(_x) << std::endl;
+ std::cerr << "(x,y) = (" << _x[0] << "," << _x[1] << ")" << std::endl;
+ }
+
+ // advanced properties
+ virtual bool constant_hessian() { return true; }
+};
+
+
+//------------------------------------------------------------------------------------------------------
+
+// Example main
+int main(void)
+{
+ std::cout << "---------- 1) Get an instance of a NProblem..." << std::endl;
+ SmallNProblem snp;
+
+ std::cout << "---------- 2) (optional for debugging) Check derivatives of problem..." << std::endl;
+ COMISO::NPDerivativeChecker npd;
+ npd.check_all(&snp);
+
+ std::cout << "---------- 3) setup list of integer variables..." << std::endl;
+ std::vector discrete_variables;
+ discrete_variables.push_back( COMISO::PairIndexVtype(0,COMISO::Integer) );
+
+ std::cout << "---------- 4) setup constraints..." << std::endl;
+ std::vector constraints;
+ // setup constraint x+y <= 6.5
+ COMISO::LinearConstraint::SVectorNC coeffs(2);
+ coeffs.coeffRef(0) = 1.0;
+ coeffs.coeffRef(1) = 1.0;
+ COMISO::LinearConstraint lc(coeffs, -6.5, COMISO::LinearConstraint::NC_LESS_EQUAL);
+ constraints.push_back(&lc);
+
+
+// check if IPOPT solver available in current configuration
+#if( COMISO_GUROBI_AVAILABLE)
+ std::cout << "---------- 5) Get GUROBI solver... " << std::endl;
+ COMISO::GUROBISolver gsol;
+
+ std::cout << "---------- 4) Solve..." << std::endl;
+
+ gsol.solve(&snp, constraints, discrete_variables);
+#endif
+
+ // check if TAO solver available in current configuration
+#if( COMISO_CPLEX_AVAILABLE)
+ std::cout << "---------- 5) Solve with CPLEX solver... " << std::endl;
+ COMISO::CPLEXSolver csol;
+
+ std::cout << "---------- 4) Solve..." << std::endl;
+
+ csol.solve(&snp, constraints, discrete_variables);
+#endif
+
+ return 0;
+}
+
diff --git a/src/external/CoMISo/Examples/small_nleast_squares/CMakeLists.txt b/src/external/CoMISo/Examples/small_nleast_squares/CMakeLists.txt
new file mode 100644
index 000000000..438cdb292
--- /dev/null
+++ b/src/external/CoMISo/Examples/small_nleast_squares/CMakeLists.txt
@@ -0,0 +1,15 @@
+include (CoMISoExample)
+
+if (WIN32)
+ acg_add_executable (small_nleast_squares WIN32 ${sources} ${headers} )
+else ()
+ acg_add_executable (small_nleast_squares ${sources} ${headers} )
+endif ()
+
+# enable rpath linking
+set_target_properties(small_nleast_squares PROPERTIES INSTALL_RPATH_USE_LINK_PATH 1)
+
+target_link_libraries (small_nleast_squares
+ CoMISo
+ ${COMISO_LINK_LIBRARIES}
+)
diff --git a/src/external/CoMISo/Examples/small_nleast_squares/main.cc b/src/external/CoMISo/Examples/small_nleast_squares/main.cc
new file mode 100644
index 000000000..4aba84765
--- /dev/null
+++ b/src/external/CoMISo/Examples/small_nleast_squares/main.cc
@@ -0,0 +1,85 @@
+/*===========================================================================*\
+ * *
+ * CoMISo *
+ * Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen *
+ * www.rwth-graphics.de *
+ * *
+ *---------------------------------------------------------------------------*
+ * This file is part of CoMISo. *
+ * *
+ * CoMISo is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * CoMISo is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with CoMISo. If not, see . *
+ * *
+\*===========================================================================*/
+
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+
+// solve least squares problem for x=1, y=2 and x-2y = 1
+
+// Example main
+int main(void)
+{
+ std::cout << "---------- 1) Get an instance of a LeastSquaresProblem..." << std::endl;
+ // number of unknowns
+ const int n = 2;
+ COMISO::LeastSquaresProblem lsqp(n);
+
+ // term0
+ COMISO::LinearConstraint::SVectorNC coeffs0(n);
+ coeffs0.coeffRef(0) = 1.0;
+ COMISO::LinearConstraint term0(coeffs0,-1.0,COMISO::NConstraintInterface::NC_EQUAL);
+ lsqp.add_term(&term0);
+
+ // term1
+ COMISO::LinearConstraint::SVectorNC coeffs1(n);
+ coeffs1.coeffRef(1) = 1.0;
+ COMISO::LinearConstraint term1(coeffs1,-2.0,COMISO::NConstraintInterface::NC_EQUAL);
+ lsqp.add_term(&term1);
+
+ // term2
+ COMISO::LinearConstraint::SVectorNC coeffs2(n);
+ coeffs2.coeffRef(0) = 1.0;
+ coeffs2.coeffRef(1) = -2.0;
+ COMISO::LinearConstraint term2(coeffs2,-1.0,COMISO::NConstraintInterface::NC_EQUAL);
+ lsqp.add_term(&term2);
+
+ std::cout << "---------- 2) (optional for debugging) Check derivatives of problem..." << std::endl;
+ COMISO::NPDerivativeChecker npd;
+ npd.check_all(&lsqp);
+
+// check if IPOPT solver available in current configuration
+#if( COMISO_IPOPT_AVAILABLE)
+ std::cout << "---------- 3) Get IPOPT solver... " << std::endl;
+ COMISO::IPOPTSolver ipsol;
+
+ std::cout << "---------- 4) Solve..." << std::endl;
+ // there are no constraints -> provide an empty vector
+ std::vector constraints;
+ ipsol.solve(&lsqp, constraints);
+#endif
+
+ std::cout << "---------- 5) Print solution..." << std::endl;
+ for( int i=0; i provide an empty vector
+ std::vector constraints;
+ ipsol.solve(&lsqp, constraints);
+#endif
+
+ std::cout << "---------- 5) Print solution..." << std::endl;
+ for( int i=0; i > A(n,n);
+ std::vector< double > x(n);
+ std::vector< double > b(n);
+
+ init_les( A, b);
+
+ // create an empty constraint matrix (will be used later)
+ gmm::row_matrix< gmm::wsvector< double > > constraints(0,n+1); //n+1 because of right hand side
+ // create an empty vector of variable indices to be rounded (will be used later)
+ std::vector< int > ids_to_round;
+
+ std::cout << A << std::endl << b << std::endl;
+
+
+ std::cout << "---------- 2) The original solution to this system is..." << std::endl;
+
+ COMISO::ConstrainedSolver cs;
+ //void solve( RMatrixT& _constraints, CMatrixT& _A, VectorT& _x, VectorT& _rhs, VectorIT& _idx_to_round, double _reg_factor = 0.0, bool _show_miso_settings = true, bool _show_timings = true );
+ //_show_miso_settings requires a QT context and hence must be false in this example
+ cs.solve( constraints, A, x, b, ids_to_round, 0.0, false, true);
+ // copy this solution for later
+ std::vector< double > org_x( x);
+ std::cout << x << std::endl;
+
+
+ std::cout << "---------- 3) Rounding: forcing the second variable to lie on an integer, changes the solution to..." << std::endl;
+ // reset system
+ init_les( A, b);
+ ids_to_round.push_back(1);
+ cs.solve( constraints, A, x, b, ids_to_round, 0.0, false, true);
+ std::cout << x << std::endl;
+
+
+ std::cout << "---------- 4) Constraining: forcing the first variable to equal the second changes the solution to..." << std::endl;
+ // reset system
+ init_les( A, b);
+ ids_to_round.clear();
+ ids_to_round.push_back(1);
+ // setup constraint x0*1+x1*0+x2*(-1)+x3*0=0
+ gmm::resize( constraints, 1, n+1);
+ constraints( 0, 0 ) = 1.0;
+ constraints( 0, 1 ) = -1.0;
+ std::cout << " the constraint equation looks like this:" << std::endl;
+ print_equations( constraints);
+ cs.solve( constraints, A, x, b, ids_to_round, 0.0, false, true);
+ std::cout << x << std::endl;
+
+ return -1;
+}
+
diff --git a/src/external/CoMISo/Examples/small_quadratic_resolve_example/CMakeLists.txt b/src/external/CoMISo/Examples/small_quadratic_resolve_example/CMakeLists.txt
new file mode 100644
index 000000000..982942ddc
--- /dev/null
+++ b/src/external/CoMISo/Examples/small_quadratic_resolve_example/CMakeLists.txt
@@ -0,0 +1,15 @@
+include (CoMISoExample)
+
+if (WIN32)
+ acg_add_executable (small_quadratic_resolve WIN32 ${sources} ${headers} )
+else ()
+ acg_add_executable (small_quadratic_resolve ${sources} ${headers} )
+endif ()
+
+# enable rpath linking
+set_target_properties(small_quadratic_resolve PROPERTIES INSTALL_RPATH_USE_LINK_PATH 1)
+
+target_link_libraries (small_quadratic_resolve
+ CoMISo
+ ${COMISO_LINK_LIBRARIES}
+)
diff --git a/src/external/CoMISo/Examples/small_quadratic_resolve_example/main.cc b/src/external/CoMISo/Examples/small_quadratic_resolve_example/main.cc
new file mode 100644
index 000000000..d0425a104
--- /dev/null
+++ b/src/external/CoMISo/Examples/small_quadratic_resolve_example/main.cc
@@ -0,0 +1,146 @@
+/*===========================================================================*\
+ * *
+ * CoMISo *
+ * Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen *
+ * www.rwth-graphics.de *
+ * *
+ *---------------------------------------------------------------------------*
+ * This file is part of CoMISo. *
+ * *
+ * CoMISo is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * CoMISo is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with CoMISo. If not, see . *
+ * *
+\*===========================================================================*/
+
+
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+
+/// function to initialize a simple system of linear equations
+template
+void init_les( MatrixT& _A, std::vector< double >& _b)
+{
+ _A(0,0) = 25 ; _A(0,1) = 0 ; _A(0,2) = -15; _A(0,3) = 0 ;
+ _A(1,0) = 0 ; _A(1,1) = 20; _A(1,2) = -8 ; _A(1,3) = -4;
+ _A(2,0) = -15 ; _A(2,1) = -8; _A(2,2) = 17 ; _A(2,3) = 0 ;
+ _A(3,0) = 0 ; _A(3,1) = -4; _A(3,2) = 0 ; _A(3,3) = 4 ;
+
+ _b[0] = 0; _b[1] = 4; _b[2] = -2; _b[3] = 0;
+}
+
+/// function to print the equations corresponding to the matrices of an equation system
+template
+void print_equations( const MatrixT& _B)
+{
+ int m = gmm::mat_nrows( _B);
+ int n = gmm::mat_ncols( _B);
+ for( int i = 0; i < m; ++i)
+ {
+ for( int j = 0; j < n-1; ++j)
+ {
+ if( _B(i,j) != 0.0)
+ std::cout << _B(i,j) << "*x" << j;
+ else
+ std::cout << " 0 ";
+ if( j < n-2 ) std::cout << " + ";
+ }
+ std::cout << " = " << -_B(i, n-1) << std::endl;
+ }
+}
+
+
+// Example main
+int main(void)
+{
+ std::cout << "---------- 1) Setup small (symmetric) test equation system Ax=b..." << std::endl;
+ int n = 4;
+ gmm::col_matrix< gmm::wsvector< double > > A(n,n);
+ std::vector< double > x(n);
+ std::vector< double > b(n);
+
+ std::vector x_bak;
+
+ std::cout << "---------- 1) Set up problem..." << std::endl;
+
+ init_les( A, b);
+
+ // create an empty constraint matrix (will be used later)
+ gmm::row_matrix< gmm::wsvector< double > > constraints(0,n+1); //n+1 because of right hand side
+ // create an empty vector of variable indices to be rounded (will be used later)
+ std::vector< int > ids_to_round;
+
+ std::cout << A << std::endl << b << std::endl;
+
+ // setup constraints
+ gmm::resize( constraints, 3, n+1);
+ constraints( 0, 0 ) = 1.0;
+ constraints( 0, 1 ) = -1.0;
+ constraints( 0, n ) = 2.0;
+ constraints( 1, 3 ) = 1.0;
+ constraints( 1, n ) = -1.0;
+ // add one redundant constraint (this will be filtered out during Gaussian elimination)
+ constraints( 2, 0 ) = 1.0;
+ constraints( 2, 1 ) = -1.0;
+ constraints( 2, n ) = 2.0;
+ std::cout << " the constraint equations looks like this:" << std::endl;
+ print_equations( constraints);
+
+ std::cout << "---------- 2) Solve full ..." << std::endl;
+ COMISO::ConstrainedSolver cs;
+ cs.solve_const( constraints, A, x, b, ids_to_round, 0.0, false, true);
+ x_bak = x;
+
+ // first test: resolve with identical rhs's
+ std::vector constraint_rhs(3);
+ std::vector b_new = b;
+ constraint_rhs[0] = -2.0;
+ constraint_rhs[1] = 1.0;
+ constraint_rhs[2] = -2.0;
+
+ std::cout << "---------- 2) Solve same rhs pre-factorized ..." << std::endl;
+ cs.resolve(x, &constraint_rhs, &b_new);
+ std::cout << "orig result: " << x_bak << std::endl;
+ std::cout << "resolve result: " << x << std::endl;
+
+ // second test: resolve with changed rhs
+ constraint_rhs[0] = 4.0;
+ constraint_rhs[1] = -2.0;
+ constraint_rhs[2] = 4.0;
+ b_new[0] = 1.0;
+ b_new[1] = -2.0;
+ b_new[2] = 3.0;
+ b_new[3] = -5.0;
+
+ std::cout << "---------- 3) Solve different rhs pre-factorized ..." << std::endl;
+ cs.resolve(x, &constraint_rhs, &b_new);
+
+
+ // solve with new factorization
+ constraints( 0, n ) = -4.0;
+ constraints( 1, n ) = 2.0;
+ constraints( 2, n ) = -4.0;
+ std::cout << "---------- 4) Solve different rhs full ..." << std::endl;
+ cs.solve_const( constraints, A, x_bak, b_new, ids_to_round, 0.0, false, true);
+
+ std::cout << "orig result (with different rhs's): " << x_bak << std::endl;
+ std::cout << "resolve result (with different rhs's): " << x << std::endl;
+
+ return 0;
+}
+
diff --git a/src/external/CoMISo/Examples/small_sparseqr/CMakeLists.txt b/src/external/CoMISo/Examples/small_sparseqr/CMakeLists.txt
new file mode 100644
index 000000000..1bc292aef
--- /dev/null
+++ b/src/external/CoMISo/Examples/small_sparseqr/CMakeLists.txt
@@ -0,0 +1,15 @@
+include (CoMISoExample)
+
+if (WIN32)
+ acg_add_executable (small_sparseqr WIN32 ${sources} ${headers} )
+else ()
+ acg_add_executable (small_sparseqr ${sources} ${headers} )
+endif ()
+
+# enable rpath linking
+set_target_properties(small_sparseqr PROPERTIES INSTALL_RPATH_USE_LINK_PATH 1)
+
+target_link_libraries (small_sparseqr
+ CoMISo
+ ${COMISO_LINK_LIBRARIES}
+)
diff --git a/src/external/CoMISo/Examples/small_sparseqr/main.cc b/src/external/CoMISo/Examples/small_sparseqr/main.cc
new file mode 100644
index 000000000..e0087abe3
--- /dev/null
+++ b/src/external/CoMISo/Examples/small_sparseqr/main.cc
@@ -0,0 +1,176 @@
+/*===========================================================================*\
+ * *
+ * CoMISo *
+ * Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen *
+ * www.rwth-graphics.de *
+ * *
+ *---------------------------------------------------------------------------*
+ * This file is part of CoMISo. *
+ * *
+ * CoMISo is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * CoMISo is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with CoMISo. If not, see . *
+ * *
+\*===========================================================================*/
+
+#include
+#include
+#include
+#include
+
+//------------------------------------------------------------------------------------------------------
+#if COMISO_SUITESPARSE_SPQR_AVAILABLE // additional spqr library required
+//------------------------------------------------------------------------------------------------------
+
+#include
+#include
+#include
+#include
+
+
+//------------------------------------------------------------------------------------------------------
+
+// Example main
+int main(void)
+{
+ std::cout << "---------- 0) Using Sparse QR for solving underdetermined equations and computing Null spaces " << std::endl;
+
+ typedef Eigen::SparseMatrix< double > SpMatrix;
+ typedef Eigen::MatrixXd DenMatrix;
+ typedef Eigen::Triplet< double > Triplet;
+
+ int dimr(4+1);
+ int dimc(4+2);
+
+ std::cout << "---------- 1) Creating matrix " << std::endl;
+ std::vector< Triplet > triplets;
+ for( int i = 0; i < dimc*dimr/2; ++i)
+ {
+ int x( rand()%(dimr-1));
+ int y( rand()%dimc);
+ double val( rand()%10);
+ //std::cerr << " setting (" << x << ", " << y << ") to " << val << std::endl;
+ triplets.push_back( Triplet( x, y, val));
+ }
+ SpMatrix A(dimr,dimc);
+ A.setFromTriplets(triplets.begin(), triplets.end());
+
+ std::cerr << DenMatrix(A) << std::endl;
+ int m = dimr;
+ int n = dimc;
+
+ if( m < n )
+ {
+ std::swap( m,n);
+ std::cerr << " ... m < n -> form transposed ..." << std::endl;
+ A = SpMatrix(A.transpose());
+ // test make also row -rank-deficinet
+ A.middleCols(n-1,1) = A.middleCols(0,1);
+ A.middleCols(0,1) = A.middleCols(n-2,1);
+ std::cerr << DenMatrix(A) << std::endl;
+ }
+
+
+ std::cerr << " ... m = " << m << "; n = " << n << std::endl;
+ std::cerr << std::endl;
+
+ std::cout << "---------- 2) Sparse QR " << std::endl;
+ COMISO::SparseQRSolver spqr;
+ SpMatrix Q,R;
+ std::vector< size_t > P;
+ int rank = spqr.factorize_system_eigen( A, Q, R, P);
+ int nullity(dimc-rank);
+ // setup permutation matrix
+ SpMatrix Pm( n, n);
+ if( !P.empty())
+ {
+ for( size_t i = 0; i < P.size(); ++i)
+ {
+ Pm.coeffRef( i, P[i]) = 1;
+ }
+ }
+
+ std::cout << "---------- 3) Result " << std::endl;
+ std::cerr << " Q " << std::endl << DenMatrix(Q) << std::endl;
+ std::cerr << " R " << std::endl << DenMatrix(R) << std::endl;
+ std::cerr << " P " << std::endl << P << std::endl;
+ std::cerr << " P matrix " << std::endl << DenMatrix(Pm) << std::endl;
+ std::cerr << " Rank " << rank << std::endl;
+ std::cerr << " Nullity " << nullity << std::endl;
+ // extract nullspace
+ SpMatrix NullSpace( Q.middleCols( std::max( 0, m-nullity), nullity));
+ std::cerr << " Nullspace " << std::endl << DenMatrix(NullSpace) << std::endl;
+ // non nullspace part of R
+ //// assuming superflous column in R is the last (if A is also row deficient)
+ //SpMatrix Rtmp(R.middleCols(0,std::min(n,n-(n-rank))).transpose());
+ //SpMatrix R1( R.transpose().middleCols(0, m-nullity));
+ SpMatrix Rtmp(R.transpose());
+ SpMatrix R1t( Rtmp.middleCols(0,m-nullity));
+ SpMatrix R1( R1t.transpose());
+ std::cerr << " Non-Nullspace R " << std::endl << DenMatrix(R1) << std::endl;
+
+
+
+ std::cout << "---------- 4) Verification " << std::endl;
+ SpMatrix reconstructedA(Q*R*Pm.transpose());
+ std::cerr << " Q orthogonal? \t " << ((fabs((Q.transpose()*Q).squaredNorm()-m) < 1e-8)?"yes":"no") << std::endl;
+ std::cerr << " A = QR? \t " << (((reconstructedA-A).squaredNorm() < 1e-8)? "yes":"no") << std::endl;
+
+
+ std::cerr << std::endl << std::endl;
+ std::cout << "---------- 5) Solving Ax=b (with x without nullspace component)" << std::endl;
+ // NOTE: A was transposed above to be m>n
+ SpMatrix b(n,1);
+ SpMatrix x(m,1);
+ for( int i = 0; i < n; ++i)
+ b.coeffRef(i,0) = rand()%10;
+ std::cerr << " ... System Ax = b .. \n";
+ std::cerr << " A " << std::endl << DenMatrix(A.transpose()) << " x " << std::endl << DenMatrix(x) << " b " << std::endl << DenMatrix(b) << std::endl;
+
+ std::cout << "---------- 5.1) test: solve using sparse QR solving .." << std::endl;
+ SpMatrix At(A.transpose());
+ spqr.solve_system_eigen( At, b, x);
+
+ std::cerr << " ... solution x .. " << std::endl;
+ std::cerr << DenMatrix(x) << std::endl;
+
+ std::cerr << " ... test: is a solution ? " << (((A.transpose()*x-b).squaredNorm()<1e-8)?"yes":"no") << std::endl;
+ std::cerr << " ... test: has nullspace component ? " << ((x.transpose()*NullSpace).squaredNorm()<1e-8?"yes":"no") << std::endl;
+ std::cerr << " ... Nullspace projections : " << (x.transpose()*NullSpace) << std::endl;
+
+ std::cout << "---------- 5.2) test: solve without nullspace .." << std::endl;
+ SpMatrix Atnull(At);
+ SpMatrix bnull(b);
+ SpMatrix xnull(m,1);
+ spqr.solve_system_eigen_min2norm( Atnull, bnull, xnull);
+ std::cerr << " ... solution x .. " << std::endl;
+ std::cerr << DenMatrix(xnull) << std::endl;
+
+ std::cerr << " ... test: is a solution ? " << (((A.transpose()*xnull-bnull).squaredNorm()<1e-8)?"yes":"no") << std::endl;
+ std::cerr << " ... test: has nullspace component ? " << ((xnull.transpose()*NullSpace).squaredNorm()<1e-8?"yes":"no") << std::endl;
+ std::cerr << " ... Nullspace projections : " << (xnull.transpose()*NullSpace) << std::endl;
+
+
+
+ return 0;
+}
+
+#else // COMISO_SUITESPARSE_SPQR_AVAILABLE
+
+int main(void)
+{
+ std::cerr << " SUITESPARSE_SPQR not available, please re-configure!\n";
+ return 0;
+}
+
+#endif // COMISO_SUITESPARSE_SPQR_AVAILABLE
+
diff --git a/src/external/CoMISo/NSolver/BoundConstraint.cc b/src/external/CoMISo/NSolver/BoundConstraint.cc
new file mode 100644
index 000000000..c0b541a8b
--- /dev/null
+++ b/src/external/CoMISo/NSolver/BoundConstraint.cc
@@ -0,0 +1,161 @@
+//=============================================================================
+//
+// CLASS BoundConstraint - IMPLEMENTATION
+//
+//=============================================================================
+
+
+//== INCLUDES =================================================================
+
+#include "BoundConstraint.hh"
+
+
+//== FORWARDDECLARATIONS ======================================================
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== CLASS DEFINITION =========================================================
+
+
+
+BoundConstraint::
+BoundConstraint(const unsigned int _var_idx, // index of variable for bound constraint
+ const double _bound, // bound: x(_var_idx) #_type, <,=,># _bound
+ const unsigned int _n, // number of unknowns in problem
+ const ConstraintType _type) // type of bound upper, lower or both (equal)
+ : NConstraintInterface(_type), idx_(_var_idx), bound_(_bound), n_(_n)
+{
+}
+
+
+//-----------------------------------------------------------------------------
+
+BoundConstraint::
+~BoundConstraint()
+{
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+int
+BoundConstraint::
+n_unknowns()
+{
+ return n_;
+}
+
+
+//-----------------------------------------------------------------------------
+
+double
+BoundConstraint::
+eval_constraint ( const double* _x )
+{
+ return _x[idx_] - bound_;
+}
+
+
+//-----------------------------------------------------------------------------
+
+void
+BoundConstraint::
+eval_gradient ( const double* _x, SVectorNC& _g )
+{
+ _g.resize(n_); _g.coeffRef(idx_) = 1.0;
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+void
+BoundConstraint::
+eval_hessian ( const double* _x, SMatrixNC& _h )
+{
+ _h.clear(); _h.resize(n_,n_);
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+bool
+BoundConstraint::
+is_linear() const
+{
+ return true;
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+bool
+BoundConstraint::
+constant_gradient() const
+{
+ return true;
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+bool
+BoundConstraint::
+constant_hessian() const
+{
+ return true;
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+unsigned int&
+BoundConstraint::
+idx()
+{
+ return idx_;
+}
+
+//-----------------------------------------------------------------------------
+
+
+double&
+BoundConstraint::
+bound()
+{
+ return bound_;
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+unsigned int&
+BoundConstraint::
+n()
+{
+ return n_;
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+void
+BoundConstraint::
+resize(const unsigned int _n)
+{
+ n_ = _n;
+}
+
+
+//=============================================================================
+} // namespace COMISO
+//=============================================================================
diff --git a/src/external/CoMISo/NSolver/BoundConstraint.hh b/src/external/CoMISo/NSolver/BoundConstraint.hh
new file mode 100644
index 000000000..c419ef141
--- /dev/null
+++ b/src/external/CoMISo/NSolver/BoundConstraint.hh
@@ -0,0 +1,84 @@
+//=============================================================================
+//
+// CLASS BoundConstraint
+//
+//=============================================================================
+
+
+#ifndef COMISO_BOUNDCONSTRAINT_HH
+#define COMISO_BOUNDCONSTRAINT_HH
+
+
+//== INCLUDES =================================================================
+
+#include
+#include "NConstraintInterface.hh"
+
+
+//== FORWARDDECLARATIONS ======================================================
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== CLASS DEFINITION =========================================================
+
+
+
+/** \class BoundConstraint
+
+ Brief Description.
+
+ A more elaborate description follows.
+*/
+class COMISODLLEXPORT BoundConstraint : public NConstraintInterface
+{
+public:
+
+// inherited from NConstraintInterface
+// typedef Eigen::SparseVector SVectorNC;
+// typedef SuperSparseMatrixT SMatrixNC;
+// // different types of constraints
+// enum ConstraintType {NC_EQUAL, NC_LESS_EQUAL, NC_GREATER_EQUAL};
+
+ /// Default constructor
+ BoundConstraint(const unsigned int _var_idx = 0, // index of variable for bound constraint
+ const double _bound = 0.0, // bound: x(_var_idx) #_type, <,=,># _bound
+ const unsigned int _n = 0, // number of unknowns in problem
+ const ConstraintType _type = NC_LESS_EQUAL); // type of bound upper, lower or both (equal)
+
+
+ /// Destructor
+ virtual ~BoundConstraint();
+
+ virtual int n_unknowns ( );
+ virtual double eval_constraint ( const double* _x );
+ virtual void eval_gradient ( const double* _x, SVectorNC& _g );
+ virtual void eval_hessian ( const double* _x, SMatrixNC& _h );
+
+ virtual bool is_linear() const;
+ virtual bool constant_gradient() const;
+ virtual bool constant_hessian () const;
+
+ // set/get values
+ unsigned int& idx();
+ double& bound();
+ unsigned int& n();
+ void resize(const unsigned int _n);
+
+private:
+ // variable idx
+ unsigned int idx_;
+ // variable bound
+ double bound_;
+ // number of unknowns
+ unsigned int n_;
+};
+
+
+//=============================================================================
+} // namespace COMISO
+//=============================================================================
+#endif // ACG_BOUNDCONSTRAINT_HH defined
+//=============================================================================
+
diff --git a/src/external/CoMISo/NSolver/COMISOSolver.cc b/src/external/CoMISo/NSolver/COMISOSolver.cc
new file mode 100644
index 000000000..ff01449a5
--- /dev/null
+++ b/src/external/CoMISo/NSolver/COMISOSolver.cc
@@ -0,0 +1,140 @@
+//=============================================================================
+//
+// CLASS COMISOSolver - IMPLEMENTATION
+//
+//=============================================================================
+
+//== INCLUDES =================================================================
+
+//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
+#include
+//=============================================================================
+
+#include
+#include "COMISOSolver.hh"
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== IMPLEMENTATION ==========================================================
+
+
+// ********** SOLVE **************** //
+void
+COMISOSolver::
+solve(NProblemInterface* _problem,
+ std::vector& _constraints,
+ std::vector& _discrete_constraints,
+ double _reg_factor,
+ bool _show_miso_settings,
+ bool _show_timings )
+{
+
+ //----------------------------------------------
+ // 1. identify integer variables
+ //----------------------------------------------
+
+ // identify integer variables
+ std::vector round_idxs;
+ for(unsigned int i=0; i<_discrete_constraints.size(); ++i)
+ switch(_discrete_constraints[i].second)
+ {
+ case Binary :
+ case Integer:
+ round_idxs.push_back(_discrete_constraints[i].first); break;
+ default : break;
+ }
+
+
+ //----------------------------------------------
+ // 2. setup constraints
+ //----------------------------------------------
+ int n = _problem->n_unknowns();
+ gmm::row_matrix< gmm::wsvector< double > > C(_constraints.size(), n+1);
+ int n_constraints = 0;
+
+ // get zero vector
+ std::vector x(n, 0.0);
+
+ for(unsigned int i=0; i<_constraints.size(); ++i)
+ if(_constraints[i]->constraint_type() == NConstraintInterface::NC_EQUAL)
+ {
+ if(!_constraints[i]->is_linear())
+ std::cerr << "Warning: COMISOSolver received a problem with non-linear constraints!!!" << std::endl;
+
+ // get linear part
+ NConstraintInterface::SVectorNC gc;
+ _constraints[i]->eval_gradient(P(x), gc);
+
+ NConstraintInterface::SVectorNC::InnerIterator v_it(gc);
+ for(; v_it; ++v_it)
+ C(n_constraints, v_it.index()) = v_it.value();
+
+ // get constant part
+ C(n_constraints, n) = _constraints[i]->eval_constraint(P(x));
+
+ // move to next constraint
+ ++n_constraints;
+ }
+
+ // resize matrix to final number of constraints
+ gmm::resize(C,n_constraints, n+1);
+
+
+ //----------------------------------------------
+ // 3. setup energy
+ //----------------------------------------------
+
+ if(!_problem->constant_hessian())
+ std::cerr << "Warning: COMISOSolver received a problem with non-constant hessian!!!" << std::endl;
+
+
+ // get hessian matrix
+ gmm::col_matrix< gmm::wsvector< double > > A(n,n);
+ NProblemInterface::SMatrixNP H;
+ _problem->eval_hessian(P(x), H);
+ for( int i=0; i rhs(_problem->n_unknowns());
+ _problem->eval_gradient(P(x), P(rhs));
+ for(unsigned int i=0; ieval_f(P(x));
+
+ //----------------------------------------------
+ // 4. solve problem
+ //----------------------------------------------
+
+ cs_.solve(C,A,x,rhs,round_idxs,
+ _reg_factor, _show_miso_settings, _show_timings);
+
+ // void solve(
+ // RMatrixT& _constraints,
+ // CMatrixT& _A,
+ // VectorT& _x,
+ // VectorT& _rhs,
+ // VectorIT& _idx_to_round,
+ // double _reg_factor = 0.0,
+ // bool _show_miso_settings = true,
+ // bool _show_timings = true );
+
+ //----------------------------------------------
+ // 5. store result
+ //----------------------------------------------
+
+ _problem->store_result(P(x));
+
+// std::cout << "COMISO Objective: " << model.get(GRB_DoubleAttr_ObjVal) << std::endl;
+}
+
+
+//=============================================================================
+} // namespace COMISO
+//=============================================================================
diff --git a/src/external/CoMISo/NSolver/COMISOSolver.hh b/src/external/CoMISo/NSolver/COMISOSolver.hh
new file mode 100644
index 000000000..f13de482b
--- /dev/null
+++ b/src/external/CoMISo/NSolver/COMISOSolver.hh
@@ -0,0 +1,85 @@
+//=============================================================================
+//
+// CLASS COMISOSolver
+//
+//=============================================================================
+
+
+#ifndef COMISO_COMISOSOLVER_HH
+#define COMISO_COMISOSOLVER_HH
+
+
+//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
+#include
+
+//== INCLUDES =================================================================
+
+#include
+#include
+#include
+#include "NProblemInterface.hh"
+#include "NConstraintInterface.hh"
+#include "VariableType.hh"
+
+
+//== FORWARDDECLARATIONS ======================================================
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== CLASS DEFINITION =========================================================
+
+
+
+/** \class NewtonSolver GUROBISolver.hh
+
+ Brief Description.
+
+ A more elaborate description follows.
+*/
+class COMISODLLEXPORT COMISOSolver
+{
+public:
+
+ typedef std::pair PairUiV;
+
+ /// Default constructor
+ COMISOSolver() {}
+
+ /// Destructor
+ ~COMISOSolver() {}
+
+ // ********** SOLVE **************** //
+ void solve(NProblemInterface* _problem, // problem instance
+ std::vector& _constraints, // linear constraints
+ std::vector& _discrete_constraints, // discrete constraint
+ double _reg_factor = 0.0, // reguluarization factor
+ bool _show_miso_settings = true, // show settings dialog
+ bool _show_timings = true ); // show timings
+
+
+ // get reference to ConstrainedSolver to manipulate parameters
+ ConstrainedSolver& solver() { return cs_;}
+
+protected:
+ double* P(std::vector& _v)
+ {
+ if( !_v.empty())
+ return ((double*)&_v[0]);
+ else
+ return 0;
+ }
+
+private:
+ ConstrainedSolver cs_;
+};
+
+
+
+//=============================================================================
+} // namespace COMISO
+//=============================================================================
+#endif // ACG_GUROBISOLVER_HH defined
+//=============================================================================
+
diff --git a/src/external/CoMISo/NSolver/CPLEXSolver.cc b/src/external/CoMISo/NSolver/CPLEXSolver.cc
new file mode 100644
index 000000000..3f51f7b08
--- /dev/null
+++ b/src/external/CoMISo/NSolver/CPLEXSolver.cc
@@ -0,0 +1,68 @@
+//=============================================================================
+//
+// CLASS CPLEXSolver - IMPLEMENTATION
+//
+//=============================================================================
+
+//== INCLUDES =================================================================
+
+//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
+#include "CPLEXSolver.hh"
+#if COMISO_CPLEX_AVAILABLE
+//=============================================================================
+
+
+#include
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== IMPLEMENTATION ==========================================================
+
+
+
+CPLEXSolver::
+CPLEXSolver()
+{
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+//void
+//CPLEXSolver::
+//set_problem_output_path( const std::string &_problem_output_path)
+//{
+// problem_output_path_ = _problem_output_path;
+//}
+//
+//
+////-----------------------------------------------------------------------------
+//
+//
+//void
+//CPLEXSolver::
+//set_problem_env_output_path( const std::string &_problem_env_output_path)
+//{
+// problem_env_output_path_ = _problem_env_output_path;
+//}
+//
+//
+////-----------------------------------------------------------------------------
+//
+//
+//void
+//CPLEXSolver::
+//set_solution_input_path(const std::string &_solution_input_path)
+//{
+// solution_input_path_ = _solution_input_path;
+//}
+
+
+//=============================================================================
+} // namespace COMISO
+//=============================================================================
+#endif // COMISO_CPLEX_AVAILABLE
+//=============================================================================
diff --git a/src/external/CoMISo/NSolver/CPLEXSolver.hh b/src/external/CoMISo/NSolver/CPLEXSolver.hh
new file mode 100644
index 000000000..9967c4def
--- /dev/null
+++ b/src/external/CoMISo/NSolver/CPLEXSolver.hh
@@ -0,0 +1,117 @@
+//=============================================================================
+//
+// CLASS CPLEXSolver
+//
+//=============================================================================
+
+
+#ifndef COMISO_CPLEXSOLVER_HH
+#define COMISO_CPLEXSOLVER_HH
+
+
+//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
+#include
+#if COMISO_CPLEX_AVAILABLE
+
+//== INCLUDES =================================================================
+
+#include
+#include
+#include
+#include "NProblemInterface.hh"
+#include "NConstraintInterface.hh"
+#include "VariableType.hh"
+
+#include
+ILOSTLBEGIN
+
+//== FORWARDDECLARATIONS ======================================================
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== CLASS DEFINITION =========================================================
+
+
+
+/** \class NewtonSolver CPLEXSolver.hh
+
+ Brief Description.
+
+ A more elaborate description follows.
+*/
+class COMISODLLEXPORT CPLEXSolver
+{
+public:
+
+ /// Default constructor
+ CPLEXSolver();
+
+ /// Destructor
+ ~CPLEXSolver() { /*env_.end();*/}
+
+ // ********** SOLVE **************** //
+ // this function has to be inline due to static linking issues
+ inline bool solve(NProblemInterface* _problem, // problem instance
+ std::vector& _constraints, // linear constraints
+ std::vector& _discrete_constraints, // discrete constraints
+ const double _time_limit = 60,
+ const bool _silent = false); // time limit in seconds
+
+ // same as above but without discrete constraints (for convenience)
+ inline bool solve(NProblemInterface* _problem, // problem instance
+ std::vector& _constraints, // linear constraints
+ const double _time_limit = 60,
+ const bool _silent = false) // time limit in seconds
+ { std::vector dc; return solve(_problem, _constraints, dc, _time_limit, _silent);}
+
+ // with handling of cone constrints
+ inline bool solve2(NProblemInterface* _problem, // problem instance
+ std::vector& _constraints, // linear constraints
+ std::vector& _discrete_constraints, // discrete constraints
+ const double _time_limit = 60,
+ const bool _silent = false); // time limit in seconds
+
+
+// void set_problem_output_path ( const std::string &_problem_output_path);
+// void set_problem_env_output_path( const std::string &_problem_env_output_path);
+// void set_solution_input_path ( const std::string &_solution_input_path);
+
+protected:
+ double* P(std::vector& _v)
+ {
+ if( !_v.empty())
+ return ((double*)&_v[0]);
+ else
+ return 0;
+ }
+
+private:
+
+ // CPLEX environment
+// IloEnv env_;
+
+ // filenames for exporting/importing gurobi solutions
+ // if string is empty nothing is imported or exported
+// std::string problem_output_path_;
+// std::string problem_env_output_path_;
+// std::string solution_input_path_;
+};
+
+
+
+//=============================================================================
+} // namespace COMISO
+
+//=============================================================================
+#endif // COMISO_CPLEX_AVAILABLE
+//=============================================================================
+#if defined(INCLUDE_TEMPLATES) && !defined(COMISO_CPLEXSOLVER_C)
+#define COMISO_CPLEXSOLVER_TEMPLATES
+#include "CPLEXSolverT.cc"
+#endif
+//=============================================================================
+#endif // ACG_CPLEXSOLVER_HH defined
+//=============================================================================
+
diff --git a/src/external/CoMISo/NSolver/CPLEXSolverT.cc b/src/external/CoMISo/NSolver/CPLEXSolverT.cc
new file mode 100644
index 000000000..adabe08e4
--- /dev/null
+++ b/src/external/CoMISo/NSolver/CPLEXSolverT.cc
@@ -0,0 +1,574 @@
+//=============================================================================
+//
+// CLASS GCPLEXSolver - IMPLEMENTATION
+//
+//=============================================================================
+
+#define COMISO_CPLEXSOLVER_C
+
+//== INCLUDES =================================================================
+
+//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
+#include "CPLEXSolver.hh"
+#include "LinearConstraint.hh"
+#include "BoundConstraint.hh"
+#include "ConeConstraint.hh"
+
+#if COMISO_CPLEX_AVAILABLE
+//=============================================================================
+
+
+#include
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== IMPLEMENTATION ==========================================================
+
+
+
+// ********** SOLVE **************** //
+bool
+CPLEXSolver::
+solve2(NProblemInterface* _problem,
+ std::vector& _constraints,
+ std::vector& _discrete_constraints,
+ const double _time_limit,
+ const bool _silent)
+{
+ try
+ {
+ //----------------------------------------------
+ // 0. set up environment
+ //----------------------------------------------
+ if(!_silent)
+ std::cerr << "cplex -> get environment...\n";
+ IloEnv env_;
+
+ if(!_silent)
+ std::cerr << "cplex -> get model...\n";
+ IloModel model(env_);
+ // model.getEnv().set(GRB_DoubleParam_TimeLimit, _time_limit);
+
+ //----------------------------------------------
+ // 1. allocate variables
+ //----------------------------------------------
+ if(!_silent)
+ std::cerr << "cplex -> allocate variables...\n";
+ // determine variable types: 0->real, 1->integer, 2->bool
+ std::vector vtypes(_problem->n_unknowns(),0);
+ for(unsigned int i=0; i<_discrete_constraints.size(); ++i)
+ if(_discrete_constraints[i].first < vtypes.size())
+ {
+ switch(_discrete_constraints[i].second)
+ {
+ case Integer: vtypes[_discrete_constraints[i].first] = 1; break;
+ case Binary : vtypes[_discrete_constraints[i].first] = 2; break;
+ default : break;
+ }
+ }
+ else
+ std::cerr << "ERROR: requested a discrete variable which is above the total number of variables"
+ << _discrete_constraints[i].first << " vs " << vtypes.size() << std::endl;
+
+ // CPLEX variables
+ std::vector vars;
+ // first all
+ for( int i=0; i<_problem->n_unknowns(); ++i)
+ switch(vtypes[i])
+ {
+ case 0 : vars.push_back( IloNumVar(env_,-IloInfinity, IloInfinity, IloNumVar::Float) ); break;
+ case 1 : vars.push_back( IloNumVar(env_, -IloIntMax, IloIntMax, IloNumVar::Int) ); break;
+ case 2 : vars.push_back( IloNumVar(env_, 0, 1, IloNumVar::Bool) ); break;
+ }
+
+
+ // Integrate new variables
+// model.update();
+
+ //----------------------------------------------
+ // 2. setup constraints
+ //----------------------------------------------
+ if(!_silent)
+ std::cerr << "cplex -> setup constraints...\n";
+
+ // get zero vector
+ std::vector x(_problem->n_unknowns(), 0.0);
+
+ for(unsigned int i=0; i<_constraints.size(); ++i)
+ {
+ if(!_constraints[i]->is_linear())
+ std::cerr << "Warning: CPLEXSolver received a problem with non-linear constraints!!!" << std::endl;
+
+ IloExpr lin_expr(env_);
+ NConstraintInterface::SVectorNC gc;
+ _constraints[i]->eval_gradient(P(x), gc);
+
+ NConstraintInterface::SVectorNC::InnerIterator v_it(gc);
+ for(; v_it; ++v_it)
+ lin_expr += vars[v_it.index()]*v_it.value();
+
+ double b = _constraints[i]->eval_constraint(P(x));
+
+
+ switch(_constraints[i]->constraint_type())
+ {
+ case NConstraintInterface::NC_EQUAL : model.add(lin_expr + b == 0); break;
+ case NConstraintInterface::NC_LESS_EQUAL : model.add(lin_expr + b <= 0); break;
+ case NConstraintInterface::NC_GREATER_EQUAL : model.add(lin_expr + b >= 0); break;
+ }
+
+ }
+
+// model.update();
+
+ //----------------------------------------------
+ // 3. setup energy
+ //----------------------------------------------
+ if(!_silent)
+ std::cerr << "cplex -> setup energy...\n";
+
+ if(!_problem->constant_hessian())
+ std::cerr << "Warning: CPLEXSolver received a problem with non-constant hessian!!!" << std::endl;
+
+// GRBQuadExpr objective;
+ IloExpr objective(env_);
+
+ // add quadratic part
+ NProblemInterface::SMatrixNP H;
+ _problem->eval_hessian(P(x), H);
+ for( int i=0; i g(_problem->n_unknowns());
+ _problem->eval_gradient(P(x), P(g));
+ for(unsigned int i=0; ieval_f(P(x));
+
+ model.add(IloMinimize(env_,objective));
+
+// model.set(GRB_IntAttr_ModelSense, 1);
+// model.setObjective(objective);
+// model.update();
+
+
+ //----------------------------------------------
+ // 4. solve problem
+ //----------------------------------------------
+ if(!_silent)
+ std::cerr << "cplex -> generate model...\n";
+ IloCplex cplex(model);
+ cplex.setParam(IloCplex::TiLim, _time_limit);
+ { // hack
+// 0 [CPX_NODESEL_DFS] Depth-first search
+// 1 [CPX_NODESEL_BESTBOUND] Best-bound search
+// 2 [CPX_NODESEL_BESTEST] Best-estimate search
+// 3 [CPX_NODESEL_BESTEST_ALT] Alternative best-estimate search
+// cplex.setParam(IloCplex::NodeSel , 0);
+ }
+ if(!_silent)
+ std::cerr << "cplex -> solve...\n";
+
+ // silent mode?
+ if(_silent)
+ cplex.setOut(env_.getNullStream());
+
+ IloBool solution_found = cplex.solve();
+
+
+// if (solution_input_path_.empty())
+// {
+// if (!problem_env_output_path_.empty())
+// {
+// std::cout << "Writing problem's environment into file \"" << problem_env_output_path_ << "\"." << std::endl;
+// model.getEnv().writeParams(problem_env_output_path_);
+// }
+// if (!problem_output_path_.empty())
+// {
+// std::cout << "Writing problem into file \"" << problem_output_path_ << "\"." << std::endl;
+// GurobiHelper::outputModelToMpsGz(model, problem_output_path_);
+// }
+//
+// model.optimize();
+// }
+// else
+// {
+// std::cout << "Reading solution from file \"" << solution_input_path_ << "\"." << std::endl;
+// }
+//
+ //----------------------------------------------
+ // 5. store result
+ //----------------------------------------------
+
+ if(solution_found != IloFalse)
+ {
+ if(!_silent)
+ std::cerr << "cplex -> store result...\n";
+ // store computed result
+ for(unsigned int i=0; istore_result(P(x));
+ }
+
+/*
+ if (solution_input_path_.empty())
+ {
+ // store computed result
+ for(unsigned int i=0; i " << oldSize << " != " << x.size() << std::endl;
+// throw std::runtime_error("Loaded solution vector doesn't have expected dimension.");
+// }
+// }
+//
+// _problem->store_result(P(x));
+//
+// // ObjVal is only available if the optimize was called.
+// if (solution_input_path_.empty())
+// std::cout << "GUROBI Objective: " << model.get(GRB_DoubleAttr_ObjVal) << std::endl;
+ return solution_found;
+ }
+ catch (IloException& e)
+ {
+ cerr << "Concert exception caught: " << e << endl;
+ return false;
+ }
+ catch (...)
+ {
+ cerr << "Unknown exception caught" << endl;
+ return false;
+ }
+
+ return false;
+}
+
+
+//-----------------------------------------------------------------------------
+
+bool
+CPLEXSolver::
+solve(NProblemInterface* _problem,
+ std::vector& _constraints,
+ std::vector& _discrete_constraints,
+ const double _time_limit,
+ const bool _silent)
+{
+ try
+ {
+ //----------------------------------------------
+ // 0. set up environment
+ //----------------------------------------------
+ if(!_silent)
+ std::cerr << "cplex -> get environment...\n";
+ IloEnv env_;
+
+ if(!_silent)
+ std::cerr << "cplex -> get model...\n";
+ IloModel model(env_);
+ // model.getEnv().set(GRB_DoubleParam_TimeLimit, _time_limit);
+
+ //----------------------------------------------
+ // 1. allocate variables and initialize limits
+ //----------------------------------------------
+ if(!_silent)
+ std::cerr << "cplex -> allocate variables...\n";
+ // determine variable types: 0->real, 1->integer, 2->bool
+ std::vector vtypes (_problem->n_unknowns(),0);
+
+ for(unsigned int i=0; i<_discrete_constraints.size(); ++i)
+ if(_discrete_constraints[i].first < vtypes.size())
+ {
+ switch(_discrete_constraints[i].second)
+ {
+ case Integer: vtypes [_discrete_constraints[i].first] = 1;
+ break;
+ case Binary : vtypes[_discrete_constraints[i].first] = 2;
+ break;
+ default : break;
+ }
+ }
+ else
+ std::cerr << "ERROR: requested a discrete variable which is above the total number of variables"
+ << _discrete_constraints[i].first << " vs " << vtypes.size() << std::endl;
+
+ // CPLEX variables
+ std::vector vars; vars.reserve(_problem->n_unknowns());
+ // first all
+ for( int i=0; i<_problem->n_unknowns(); ++i)
+ switch(vtypes[i])
+ {
+ case 0 : vars.push_back( IloNumVar(env_,-IloInfinity, IloInfinity, IloNumVar::Float) ); break;
+ case 1 : vars.push_back( IloNumVar(env_, -IloIntMax, IloIntMax, IloNumVar::Int) ); break;
+ case 2 : vars.push_back( IloNumVar(env_, 0, 1, IloNumVar::Bool) ); break;
+ }
+
+
+
+ // Integrate new variables
+// model.update();
+
+ //----------------------------------------------
+ // 2. setup constraints
+ //----------------------------------------------
+ if(!_silent)
+ std::cerr << "cplex -> setup constraints...\n";
+
+ // get zero vector
+ std::vector x(_problem->n_unknowns(), 0.0);
+
+ // handle constraints depending on their tyep
+ for(unsigned int i=0; i<_constraints.size(); ++i)
+ {
+ // is linear constraint?
+ LinearConstraint* lin_ptr = dynamic_cast(_constraints[i]);
+ if(lin_ptr)
+ {
+ IloExpr lin_expr(env_);
+ NConstraintInterface::SVectorNC gc;
+ _constraints[i]->eval_gradient(P(x), gc);
+
+ NConstraintInterface::SVectorNC::InnerIterator v_it(gc);
+ for(; v_it; ++v_it)
+ lin_expr += vars[v_it.index()]*v_it.value();
+
+ double b = _constraints[i]->eval_constraint(P(x));
+
+
+ switch(_constraints[i]->constraint_type())
+ {
+ case NConstraintInterface::NC_EQUAL : model.add(lin_expr + b == 0); break;
+ case NConstraintInterface::NC_LESS_EQUAL : model.add(lin_expr + b <= 0); break;
+ case NConstraintInterface::NC_GREATER_EQUAL : model.add(lin_expr + b >= 0); break;
+ }
+ }
+ else
+ {
+ BoundConstraint* bnd_ptr = dynamic_cast(_constraints[i]);
+ if(bnd_ptr)
+ {
+ switch(bnd_ptr->constraint_type())
+ {
+ case NConstraintInterface::NC_EQUAL : vars[bnd_ptr->idx()].setBounds(bnd_ptr->bound(), bnd_ptr->bound()); break;
+ case NConstraintInterface::NC_LESS_EQUAL : vars[bnd_ptr->idx()].setUB(bnd_ptr->bound()); break;
+ case NConstraintInterface::NC_GREATER_EQUAL : vars[bnd_ptr->idx()].setLB(bnd_ptr->bound()); break;
+ }
+ }
+ else
+ {
+ ConeConstraint* cone_ptr = dynamic_cast(_constraints[i]);
+ if(cone_ptr)
+ {
+ IloExpr soc_lhs(env_);
+ IloExpr soc_rhs(env_);
+
+ soc_rhs= 0.5*cone_ptr->c()*vars[cone_ptr->i()]*vars[cone_ptr->i()];
+
+ NConstraintInterface::SMatrixNC::iterator q_it = cone_ptr->Q().begin();
+ for(; q_it != cone_ptr->Q().end(); ++q_it)
+ {
+ soc_lhs += 0.5*(*q_it)*vars[q_it.col()]*vars[q_it.row()];
+ }
+
+ model.add(soc_lhs <= soc_rhs);
+ }
+ else
+ std::cerr << "Warning: CPLEXSolver received a constraint of unknow type!!! -> skipping it" << std::endl;
+ }
+ }
+ }
+// model.update();
+
+ //----------------------------------------------
+ // 3. setup energy
+ //----------------------------------------------
+ if(!_silent)
+ std::cerr << "cplex -> setup energy...\n";
+
+ if(!_problem->constant_hessian())
+ std::cerr << "Warning: CPLEXSolver received a problem with non-constant hessian!!!" << std::endl;
+
+// GRBQuadExpr objective;
+ IloExpr objective(env_);
+
+ // add quadratic part
+ NProblemInterface::SMatrixNP H;
+ _problem->eval_hessian(P(x), H);
+ for( int i=0; i g(_problem->n_unknowns());
+ _problem->eval_gradient(P(x), P(g));
+ for(unsigned int i=0; ieval_f(P(x));
+
+ model.add(IloMinimize(env_,objective));
+
+// model.set(GRB_IntAttr_ModelSense, 1);
+// model.setObjective(objective);
+// model.update();
+
+
+ //----------------------------------------------
+ // 4. solve problem
+ //----------------------------------------------
+ if(!_silent)
+ std::cerr << "cplex -> generate model...\n";
+ IloCplex cplex(model);
+ cplex.setParam(IloCplex::TiLim, _time_limit);
+ { // hack
+// 0 [CPX_NODESEL_DFS] Depth-first search
+// 1 [CPX_NODESEL_BESTBOUND] Best-bound search
+// 2 [CPX_NODESEL_BESTEST] Best-estimate search
+// 3 [CPX_NODESEL_BESTEST_ALT] Alternative best-estimate search
+// cplex.setParam(IloCplex::NodeSel , 0);
+ }
+ if(!_silent)
+ std::cerr << "cplex -> solve...\n";
+
+ // silent mode?
+ if(_silent)
+ cplex.setOut(env_.getNullStream());
+
+ IloBool solution_found = cplex.solve();
+
+
+// if (solution_input_path_.empty())
+// {
+// if (!problem_env_output_path_.empty())
+// {
+// std::cout << "Writing problem's environment into file \"" << problem_env_output_path_ << "\"." << std::endl;
+// model.getEnv().writeParams(problem_env_output_path_);
+// }
+// if (!problem_output_path_.empty())
+// {
+// std::cout << "Writing problem into file \"" << problem_output_path_ << "\"." << std::endl;
+// GurobiHelper::outputModelToMpsGz(model, problem_output_path_);
+// }
+//
+// model.optimize();
+// }
+// else
+// {
+// std::cout << "Reading solution from file \"" << solution_input_path_ << "\"." << std::endl;
+// }
+//
+ //----------------------------------------------
+ // 5. store result
+ //----------------------------------------------
+
+ if(solution_found != IloFalse)
+ {
+ if(!_silent)
+ std::cerr << "cplex -> store result...\n";
+ // store computed result
+ for(unsigned int i=0; istore_result(P(x));
+ }
+
+/*
+ if (solution_input_path_.empty())
+ {
+ // store computed result
+ for(unsigned int i=0; i " << oldSize << " != " << x.size() << std::endl;
+// throw std::runtime_error("Loaded solution vector doesn't have expected dimension.");
+// }
+// }
+//
+// _problem->store_result(P(x));
+//
+// // ObjVal is only available if the optimize was called.
+// if (solution_input_path_.empty())
+// std::cout << "GUROBI Objective: " << model.get(GRB_DoubleAttr_ObjVal) << std::endl;
+ return solution_found;
+ }
+ catch (IloException& e)
+ {
+ cerr << "Concert exception caught: " << e << endl;
+ return false;
+ }
+ catch (...)
+ {
+ cerr << "Unknown exception caught" << endl;
+ return false;
+ }
+
+ return false;
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+//void
+//CPLEXSolver::
+//set_problem_output_path( const std::string &_problem_output_path)
+//{
+// problem_output_path_ = _problem_output_path;
+//}
+//
+//
+////-----------------------------------------------------------------------------
+//
+//
+//void
+//CPLEXSolver::
+//set_problem_env_output_path( const std::string &_problem_env_output_path)
+//{
+// problem_env_output_path_ = _problem_env_output_path;
+//}
+//
+//
+////-----------------------------------------------------------------------------
+//
+//
+//void
+//CPLEXSolver::
+//set_solution_input_path(const std::string &_solution_input_path)
+//{
+// solution_input_path_ = _solution_input_path;
+//}
+
+
+//=============================================================================
+} // namespace COMISO
+//=============================================================================
+#endif // COMISO_CPLEX_AVAILABLE
+//=============================================================================
diff --git a/src/external/CoMISo/NSolver/ConeConstraint.cc b/src/external/CoMISo/NSolver/ConeConstraint.cc
new file mode 100644
index 000000000..25b7fecc2
--- /dev/null
+++ b/src/external/CoMISo/NSolver/ConeConstraint.cc
@@ -0,0 +1,134 @@
+//=============================================================================
+//
+// CLASS ConerConstraint
+//
+//=============================================================================
+
+
+#ifndef COMISO_CONECONSTRAINT_CC
+#define COMISO_CONECONSTRAINT_CC
+
+
+//== INCLUDES =================================================================
+
+#include
+#include "NConstraintInterface.hh"
+#include "ConeConstraint.hh"
+
+//== FORWARDDECLARATIONS ======================================================
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== Implementation =========================================================
+
+
+/// Default constructor
+ConeConstraint::ConeConstraint()
+: NConstraintInterface(NConstraintInterface::NC_GREATER_EQUAL)
+{
+ Q_.clear();
+ i_ = 1.0;
+ c_ = 1.0;
+}
+
+// cone constraint of the form -> 0.5*(c_ * x(i_)^2 - x^T Q_ x) >= 0
+ConeConstraint::ConeConstraint(const double _c, const int _i, const SMatrixNC& _Q)
+: NConstraintInterface(NConstraintInterface::NC_GREATER_EQUAL),
+ c_(_c), i_(_i), Q_(_Q)
+{
+}
+
+/// Destructor
+ConeConstraint::~ConeConstraint() {}
+
+int ConeConstraint::n_unknowns()
+{
+ return Q_.cols();
+}
+
+void ConeConstraint::resize(const unsigned int _n)
+{
+ Q_.resize(_n,_n);
+}
+
+void ConeConstraint::clear()
+{
+ Q_.clear();
+}
+
+
+const ConeConstraint::SMatrixNC& ConeConstraint::Q() const
+{
+ return Q_;
+}
+ConeConstraint::SMatrixNC& ConeConstraint::Q()
+{
+ return Q_;
+}
+
+const int& ConeConstraint::i() const
+{
+ return i_;
+}
+int& ConeConstraint::i()
+{
+ return i_;
+}
+
+const double& ConeConstraint::c() const
+{
+ return c_;
+}
+double& ConeConstraint::c()
+{
+ return c_;
+}
+
+double ConeConstraint::eval_constraint ( const double* _x )
+{
+ // cone constraint of the form -> 0.5*(c_ * x(i_)^2 - x^T Q_ x) >= 0
+ double v = c_*_x[i_]*_x[i_];
+
+ SMatrixNC::iterator m_it = Q_.begin();
+ SMatrixNC::iterator m_end = Q_.end();
+
+ for(; m_it != m_end; ++m_it)
+ {
+ v -= (*m_it)*_x[m_it.row()]*_x[m_it.col()];
+ }
+
+ return 0.5*v;
+}
+
+void ConeConstraint::eval_gradient( const double* _x, SVectorNC& _g )
+{
+ _g.setZero();
+ _g.resize(Q_.rows());
+
+ SMatrixNC::iterator m_it = Q_.begin();
+ SMatrixNC::iterator m_end = Q_.end();
+
+ for(; m_it != m_end; ++m_it)
+ {
+ _g.coeffRef(m_it.row()) -= (*m_it)*_x[m_it.col()];
+ }
+
+ _g.coeffRef(i_) += c_*_x[i_];
+}
+
+void ConeConstraint::eval_hessian ( const double* _x, SMatrixNC& _h )
+{
+ _h = Q_;
+ _h.scale(-1.0);
+ _h(i_,i_) += c_;
+}
+
+
+//=============================================================================
+} // namespace COMISO
+//=============================================================================
+#endif // ACG_ConeConstraint_HH defined
+//=============================================================================
+
diff --git a/src/external/CoMISo/NSolver/ConeConstraint.hh b/src/external/CoMISo/NSolver/ConeConstraint.hh
new file mode 100644
index 000000000..3f556d310
--- /dev/null
+++ b/src/external/CoMISo/NSolver/ConeConstraint.hh
@@ -0,0 +1,90 @@
+//=============================================================================
+//
+// CLASS ConeConstraint
+//
+//=============================================================================
+
+
+#ifndef COMISO_CONECONSTRAINT_HH
+#define COMISO_CONECONSTRAINT_HH
+
+
+//== INCLUDES =================================================================
+
+#include
+#include "NConstraintInterface.hh"
+#include
+
+
+//== FORWARDDECLARATIONS ======================================================
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== CLASS DEFINITION =========================================================
+
+
+class COMISODLLEXPORT ConeConstraint : public NConstraintInterface
+{
+public:
+
+ // sparse vector type
+ typedef NConstraintInterface::SVectorNC SVectorNC;
+ typedef NConstraintInterface::SMatrixNC SMatrixNC;
+
+ /// Default constructor
+ ConeConstraint();
+
+ // cone constraint of the form -> 0.5*(c_ * x(i_)^2 - x^T Q_ x) >= 0
+ ConeConstraint(const double _c, const int _i, const SMatrixNC& _Q);
+
+ /// Destructor
+ virtual ~ConeConstraint();
+
+ virtual int n_unknowns();
+
+ // resize coefficient vector = #unknowns
+ void resize(const unsigned int _n);
+
+ // clear to zero constraint 0 =_type 0
+ void clear();
+
+ const double& c() const;
+ double& c();
+
+ const int& i() const;
+ int& i();
+
+ const SMatrixNC& Q() const;
+ SMatrixNC& Q();
+
+
+ virtual double eval_constraint ( const double* _x );
+
+ virtual void eval_gradient( const double* _x, SVectorNC& _g );
+
+ virtual void eval_hessian ( const double* _x, SMatrixNC& _h );
+
+ virtual bool is_linear() const { return false;}
+ virtual bool constant_gradient() const { return false;}
+ virtual bool constant_hessian () const { return true;}
+
+private:
+
+ // cone constraint of the form -> 0.5*(c_ * x(i_)^2 - x^T Q_ x) >= 0
+ double c_;
+ int i_;
+ SMatrixNC Q_;
+};
+
+
+//=============================================================================
+} // namespace COMISO
+//=============================================================================
+// support std vectors
+EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(COMISO::ConeConstraint);
+//=============================================================================
+#endif // ACG_CONECONSTRAINT_HH defined
+//=============================================================================
+
diff --git a/src/external/CoMISo/NSolver/GUROBISolver.cc b/src/external/CoMISo/NSolver/GUROBISolver.cc
new file mode 100644
index 000000000..313b44193
--- /dev/null
+++ b/src/external/CoMISo/NSolver/GUROBISolver.cc
@@ -0,0 +1,254 @@
+//=============================================================================
+//
+// CLASS GUROBISolver - IMPLEMENTATION
+//
+//=============================================================================
+
+//== INCLUDES =================================================================
+
+//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
+#include
+#if COMISO_GUROBI_AVAILABLE
+//=============================================================================
+
+
+#include "GUROBISolver.hh"
+
+#include
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== IMPLEMENTATION ==========================================================
+
+
+
+GUROBISolver::
+GUROBISolver()
+{
+}
+
+//-----------------------------------------------------------------------------
+
+// ********** SOLVE **************** //
+bool
+GUROBISolver::
+solve(NProblemInterface* _problem,
+ std::vector& _constraints,
+ std::vector& _discrete_constraints,
+ const double _time_limit)
+{
+ try
+ {
+ //----------------------------------------------
+ // 0. set up environment
+ //----------------------------------------------
+
+ GRBEnv env = GRBEnv();
+ GRBModel model = GRBModel(env);
+
+ model.getEnv().set(GRB_DoubleParam_TimeLimit, _time_limit);
+
+
+ //----------------------------------------------
+ // 1. allocate variables
+ //----------------------------------------------
+
+ // determine variable types: 0->real, 1->integer, 2->bool
+ std::vector vtypes(_problem->n_unknowns(),0);
+ for(unsigned int i=0; i<_discrete_constraints.size(); ++i)
+ switch(_discrete_constraints[i].second)
+ {
+ case Integer: vtypes[_discrete_constraints[i].first] = 1; break;
+ case Binary : vtypes[_discrete_constraints[i].first] = 2; break;
+ default : break;
+ }
+
+ // GUROBI variables
+ std::vector vars;
+ // first all
+ for( int i=0; i<_problem->n_unknowns(); ++i)
+ switch(vtypes[i])
+ {
+ case 0 : vars.push_back( model.addVar(-GRB_INFINITY, GRB_INFINITY, 0.0, GRB_CONTINUOUS) ); break;
+ case 1 : vars.push_back( model.addVar(-GRB_INFINITY, GRB_INFINITY, 0.0, GRB_INTEGER ) ); break;
+ case 2 : vars.push_back( model.addVar(-GRB_INFINITY, GRB_INFINITY, 0.0, GRB_BINARY ) ); break;
+ }
+
+
+ // Integrate new variables
+ model.update();
+
+ //----------------------------------------------
+ // 2. setup constraints
+ //----------------------------------------------
+
+ // get zero vector
+ std::vector x(_problem->n_unknowns(), 0.0);
+
+ for(unsigned int i=0; i<_constraints.size(); ++i)
+ {
+ if(!_constraints[i]->is_linear())
+ std::cerr << "Warning: GUROBISolver received a problem with non-linear constraints!!!" << std::endl;
+
+ GRBLinExpr lin_expr;
+ NConstraintInterface::SVectorNC gc;
+ _constraints[i]->eval_gradient(P(x), gc);
+
+ NConstraintInterface::SVectorNC::InnerIterator v_it(gc);
+ for(; v_it; ++v_it)
+// lin_expr += v_it.value()*vars[v_it.index()];
+ lin_expr = lin_expr + vars[v_it.index()]*v_it.value();
+
+ double b = _constraints[i]->eval_constraint(P(x));
+
+ switch(_constraints[i]->constraint_type())
+ {
+ case NConstraintInterface::NC_EQUAL : model.addConstr(lin_expr + b == 0); break;
+ case NConstraintInterface::NC_LESS_EQUAL : model.addConstr(lin_expr + b <= 0); break;
+ case NConstraintInterface::NC_GREATER_EQUAL : model.addConstr(lin_expr + b >= 0); break;
+ }
+ }
+ model.update();
+
+ //----------------------------------------------
+ // 3. setup energy
+ //----------------------------------------------
+
+ if(!_problem->constant_hessian())
+ std::cerr << "Warning: GUROBISolver received a problem with non-constant hessian!!!" << std::endl;
+
+ GRBQuadExpr objective;
+
+ // add quadratic part
+ NProblemInterface::SMatrixNP H;
+ _problem->eval_hessian(P(x), H);
+ for( int i=0; i g(_problem->n_unknowns());
+ _problem->eval_gradient(P(x), P(g));
+ for(unsigned int i=0; ieval_f(P(x));
+
+ model.set(GRB_IntAttr_ModelSense, 1);
+ model.setObjective(objective);
+ model.update();
+
+
+ //----------------------------------------------
+ // 4. solve problem
+ //----------------------------------------------
+
+
+ if (solution_input_path_.empty())
+ {
+ if (!problem_env_output_path_.empty())
+ {
+ std::cout << "Writing problem's environment into file \"" << problem_env_output_path_ << "\"." << std::endl;
+ model.getEnv().writeParams(problem_env_output_path_);
+ }
+ if (!problem_output_path_.empty())
+ {
+ std::cout << "Writing problem into file \"" << problem_output_path_ << "\"." << std::endl;
+ GurobiHelper::outputModelToMpsGz(model, problem_output_path_);
+ }
+
+ model.optimize();
+ }
+ else
+ {
+ std::cout << "Reading solution from file \"" << solution_input_path_ << "\"." << std::endl;
+ }
+
+ //----------------------------------------------
+ // 5. store result
+ //----------------------------------------------
+
+ if (solution_input_path_.empty())
+ {
+ // store computed result
+ for(unsigned int i=0; i " << oldSize << " != " << x.size() << std::endl;
+ throw std::runtime_error("Loaded solution vector doesn't have expected dimension.");
+ }
+ }
+
+ _problem->store_result(P(x));
+
+ // ObjVal is only available if the optimize was called.
+ if (solution_input_path_.empty())
+ std::cout << "GUROBI Objective: " << model.get(GRB_DoubleAttr_ObjVal) << std::endl;
+ return true;
+ }
+ catch(GRBException& e)
+ {
+ std::cout << "Error code = " << e.getErrorCode() << std::endl;
+ std::cout << e.getMessage() << std::endl;
+ return false;
+ }
+ catch(...)
+ {
+ std::cout << "Exception during optimization" << std::endl;
+ return false;
+ }
+
+ return false;
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+void
+GUROBISolver::
+set_problem_output_path( const std::string &_problem_output_path)
+{
+ problem_output_path_ = _problem_output_path;
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+void
+GUROBISolver::
+set_problem_env_output_path( const std::string &_problem_env_output_path)
+{
+ problem_env_output_path_ = _problem_env_output_path;
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+void
+GUROBISolver::
+set_solution_input_path(const std::string &_solution_input_path)
+{
+ solution_input_path_ = _solution_input_path;
+}
+
+
+//=============================================================================
+} // namespace COMISO
+//=============================================================================
+#endif // COMISO_GUROBI_AVAILABLE
+//=============================================================================
diff --git a/src/external/CoMISo/NSolver/GUROBISolver.hh b/src/external/CoMISo/NSolver/GUROBISolver.hh
new file mode 100644
index 000000000..caf2c7a26
--- /dev/null
+++ b/src/external/CoMISo/NSolver/GUROBISolver.hh
@@ -0,0 +1,92 @@
+//=============================================================================
+//
+// CLASS GUROBISolver
+//
+//=============================================================================
+
+
+#ifndef COMISO_GUROBISOLVER_HH
+#define COMISO_GUROBISOLVER_HH
+
+
+//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
+#include
+#if COMISO_GUROBI_AVAILABLE
+
+//== INCLUDES =================================================================
+
+#include
+#include
+#include
+#include "NProblemInterface.hh"
+#include "NConstraintInterface.hh"
+#include "VariableType.hh"
+#include "GurobiHelper.hh"
+
+#include
+
+//== FORWARDDECLARATIONS ======================================================
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== CLASS DEFINITION =========================================================
+
+
+
+/** \class NewtonSolver GUROBISolver.hh
+
+ Brief Description.
+
+ A more elaborate description follows.
+*/
+class COMISODLLEXPORT GUROBISolver
+{
+public:
+
+ /// Default constructor
+ GUROBISolver();
+
+ /// Destructor
+ ~GUROBISolver() {}
+
+ // ********** SOLVE **************** //
+ bool solve(NProblemInterface* _problem, // problem instance
+ std::vector& _constraints, // linear constraints
+ std::vector& _discrete_constraints, // discrete constraints
+ const double _time_limit = 60 ); // time limit in seconds
+
+ void set_problem_output_path ( const std::string &_problem_output_path);
+ void set_problem_env_output_path( const std::string &_problem_env_output_path);
+ void set_solution_input_path ( const std::string &_solution_input_path);
+
+protected:
+ double* P(std::vector& _v)
+ {
+ if( !_v.empty())
+ return ((double*)&_v[0]);
+ else
+ return 0;
+ }
+
+private:
+
+ // filenames for exporting/importing gurobi solutions
+ // if string is empty nothing is imported or exported
+ std::string problem_output_path_;
+ std::string problem_env_output_path_;
+ std::string solution_input_path_;
+};
+
+
+
+//=============================================================================
+} // namespace COMISO
+
+//=============================================================================
+#endif // COMISO_GUROBI_AVAILABLE
+//=============================================================================
+#endif // ACG_GUROBISOLVER_HH defined
+//=============================================================================
+
diff --git a/src/external/CoMISo/NSolver/GurobiHelper.cc b/src/external/CoMISo/NSolver/GurobiHelper.cc
new file mode 100644
index 000000000..0bfa41134
--- /dev/null
+++ b/src/external/CoMISo/NSolver/GurobiHelper.cc
@@ -0,0 +1,175 @@
+/*
+ * GurobiHelper.cc
+ *
+ * Created on: Jan 4, 2012
+ * Author: ebke
+ */
+
+#include "GurobiHelper.hh"
+
+#if (COMISO_GUROBI_AVAILABLE && COMISO_BOOST_AVAILABLE)
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+
+#define OUTPUT_UNCOMPRESSED_WITH_CONSTANT_COMMENT 0
+#define OUTPUT_CONSTANT_AS_CONT 1
+
+namespace COMISO {
+
+/**
+ * Helper class that ensures exception safe deletion
+ * of a temporary file.
+ */
+class TempFileGuard {
+ public:
+ TempFileGuard(const std::string &_filePath) : filePath_(_filePath) {
+ }
+
+ ~TempFileGuard() {
+ if (boost::filesystem::exists(filePath_))
+ boost::filesystem::remove(filePath_);
+ }
+
+ const boost::filesystem::path &filePath() const { return filePath_; };
+
+ private:
+ boost::filesystem::path filePath_;
+};
+
+static void moveConstantTermIntoConstrainedVariable(GRBModel &model) {
+ const double constantTerm = model.getObjective().getLinExpr().getConstant();
+ //tmpModel.getObjective().addConstant(-constantTerm);
+ model.getObjective() -= constantTerm;
+#if OUTPUT_CONSTANT_AS_CONT
+ model.addVar(constantTerm, constantTerm, 1, GRB_CONTINUOUS, "MIQ_synthetic_constant");
+#else
+ model.addVar(1, 1, constantTerm, GRB_INTEGER, "MIQ_synthetic_constant");
+#endif
+}
+
+static void copyFile(const char *from, const char *to) {
+ FILE *inF = fopen(from, "r");
+ FILE *outF = fopen(to, "w");
+
+ const int bufsize = 4096;
+ unsigned char buffer[bufsize];
+
+ do {
+ size_t readBytes = fread(buffer, 1, bufsize, inF);
+ fwrite(buffer, 1, readBytes, outF);
+ } while(!feof(inF));
+
+ fclose(inF);
+ fclose(outF);
+}
+
+/**
+ * WARNING: Never call outputModelToMpsGz and importInitialSolutionIntoModel
+ * on the same model. Both try to move the constant term into a variable and
+ * consequently, the second attempt to do so will fail.
+ */
+void GurobiHelper::outputModelToMpsGz(GRBModel &model, const std::string &problem_output_path_) {
+#if OUTPUT_UNCOMPRESSED_WITH_CONSTANT_COMMENT
+ boost::scoped_ptr tempFileGuard;
+ {
+ QTemporaryFile tempFile("XXXXXX.mps");
+ tempFile.setAutoRemove(false);
+ tempFile.open();
+
+ // In order to minimize the likelihood of race conditions,
+ // we initialize tempFileGuard right here.
+ tempFileGuard.reset(new TempFileGuard(QFileInfo(tempFile).absoluteFilePath().toStdString()));
+ tempFile.close();
+ }
+
+ const std::string fileName = tempFileGuard->filePath().string();
+
+ model.write(fileName);
+ const double constantTerm = model.getObjective().getLinExpr().getConstant();
+
+ FILE *inF = fopen(fileName.c_str(), "r");
+ FILE *outF = fopen(problem_output_path_.c_str(), "w");
+
+ fprintf(outF, "* Constant Term: %.16e\n", constantTerm);
+ const int bufsize = 4096;
+ unsigned char buffer[bufsize];
+ int readBytes;
+ do {
+ readBytes = fread(buffer, 1, bufsize, inF);
+ fwrite(buffer, 1, readBytes, outF);
+ } while(!feof(inF));
+ fclose(inF);
+ fclose(outF);
+#else
+ GRBModel tmpModel(model);
+
+ moveConstantTermIntoConstrainedVariable(tmpModel);
+
+ tmpModel.update();
+ tmpModel.write(problem_output_path_);
+#endif
+}
+
+/**
+ * WARNING: Never call outputModelToMpsGz and importInitialSolutionIntoModel
+ * on the same model. Both try to move the constant term into a variable and
+ * consequently, the second attempt to do so will fail.
+ */
+void GurobiHelper::importInitialSolutionIntoModel(GRBModel &model, const std::string &solution_path_) {
+ boost::scoped_ptr tempFileGuard;
+ {
+ QTemporaryFile tempFile("XXXXXX.mst");
+ tempFile.setAutoRemove(false);
+ tempFile.open();
+
+ // In order to minimize the likelihood of race conditions,
+ // we initialize tempFileGuard right here.
+ tempFileGuard.reset(new TempFileGuard(QFileInfo(tempFile).absoluteFilePath().toStdString()));
+ tempFile.close();
+ }
+
+ const std::string fileName = tempFileGuard->filePath().string();
+
+ copyFile(solution_path_.c_str(), fileName.c_str());
+
+ //moveConstantTermIntoConstrainedVariable(model);
+ const double constantTerm = model.getObjective().getLinExpr().getConstant();
+ model.addVar(constantTerm, constantTerm, 0, GRB_CONTINUOUS, "MIQ_synthetic_constant");
+
+ model.update();
+ model.read(fileName);
+ model.update();
+}
+
+void GurobiHelper::readSolutionVectorFromSOL(std::vector &out_solution_, const std::string &fileName_) {
+ std::ifstream solFile(fileName_.c_str());
+ //if (!solFile.good())
+ // throw std::runtime_error("Unable to open file \"" + fileName + "\".");
+
+ static const boost::regex commentRe("\\s*#", boost::regex_constants::perl);
+ static const boost::regex variableRe("\\s*(\\S+)\\s+([-+]?[0-9]*\\.?[0-9]+(?:[eE][-+]?[0-9]+)?)", boost::regex_constants::perl);
+
+ std::string line;
+ while (solFile) {
+ std::getline(solFile, line);
+ if (boost::regex_search(line, commentRe, boost::match_continuous)) continue;
+ boost::smatch match;
+ if (boost::regex_search(line, match, variableRe, boost::match_continuous) && match[1] != "MIQ_synthetic_constant") {
+ out_solution_.push_back(boost::lexical_cast(match[2]));
+ }
+ }
+}
+
+} /* namespace COMISO */
+
+#endif
diff --git a/src/external/CoMISo/NSolver/GurobiHelper.hh b/src/external/CoMISo/NSolver/GurobiHelper.hh
new file mode 100644
index 000000000..58d872b8b
--- /dev/null
+++ b/src/external/CoMISo/NSolver/GurobiHelper.hh
@@ -0,0 +1,51 @@
+/*
+ * GurobiHelper.hh
+ *
+ * Created on: Jan 4, 2012
+ * Author: ebke
+ */
+
+
+#ifndef GUROBIHELPER_HH_
+#define GUROBIHELPER_HH_
+
+//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
+#include
+#if (COMISO_GUROBI_AVAILABLE && COMISO_BOOST_AVAILABLE)
+//=============================================================================
+
+#include
+#include
+#include
+
+namespace COMISO {
+
+class GurobiHelper {
+ public:
+
+ /**
+ * WARNING: Never call outputModelToMpsGz and importInitialSolutionIntoModel
+ * on the same model. Both try to move the constant term into a variable and
+ * consequently, the second attempt to do so will fail.
+ */
+ static void outputModelToMpsGz(GRBModel &model, const std::string &problem_output_path_);
+
+ /**
+ * WARNING: Never call outputModelToMpsGz and importInitialSolutionIntoModel
+ * on the same model. Both try to move the constant term into a variable and
+ * consequently, the second attempt to do so will fail.
+ */
+ static void importInitialSolutionIntoModel(GRBModel &model, const std::string &solution_path_);
+
+ /**
+ * Reads the solution vector from a SOL file and appends it to
+ * out_solution_.
+ */
+ static void readSolutionVectorFromSOL(std::vector &out_solution_, const std::string &fileName_);
+};
+
+} /* namespace COMISO */
+#endif /* COMISO_GUROBI_AVAILABLE */
+#endif /* GUROBIHELPER_HH_ */
+
+
diff --git a/src/external/CoMISo/NSolver/IPOPTSolver.cc b/src/external/CoMISo/NSolver/IPOPTSolver.cc
new file mode 100644
index 000000000..8fd273dea
--- /dev/null
+++ b/src/external/CoMISo/NSolver/IPOPTSolver.cc
@@ -0,0 +1,1153 @@
+//=============================================================================
+//
+// CLASS IPOPTSolver - IMPLEMENTATION
+//
+//=============================================================================
+
+//== INCLUDES =================================================================
+
+//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
+#include
+#if COMISO_IPOPT_AVAILABLE
+//=============================================================================
+
+
+#include "IPOPTSolver.hh"
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== IMPLEMENTATION ==========================================================
+
+
+// Constructor
+IPOPTSolver::
+IPOPTSolver()
+{
+ // Create an instance of the IpoptApplication
+ app_ = IpoptApplicationFactory();
+
+ // Switch to HSL if available in Comiso
+ #if COMISO_HSL_AVAILABLE
+ app_->Options()->SetStringValue("linear_solver", "ma57");
+ #endif
+
+ // Restrict memory to be able to run larger problems on windows
+ // with the default mumps solver
+ #ifdef WIN32
+ app_->Options()->SetIntegerValue("mumps_mem_percent", 5);
+ #endif
+
+ // set default parameters
+ app_->Options()->SetIntegerValue("max_iter", 100);
+ // app->Options()->SetStringValue("derivative_test", "second-order");
+ // app->Options()->SetIntegerValue("print_level", 0);
+ // app->Options()->SetStringValue("expect_infeasible_problem", "yes");
+
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+
+int
+IPOPTSolver::
+solve(NProblemInterface* _problem, const std::vector& _constraints)
+{
+ //----------------------------------------------------------------------------
+ // 1. Create an instance of IPOPT NLP
+ //----------------------------------------------------------------------------
+ Ipopt::SmartPtr np = new NProblemIPOPT(_problem, _constraints);
+ NProblemIPOPT* np2 = dynamic_cast (Ipopt::GetRawPtr(np));
+
+ //----------------------------------------------------------------------------
+ // 2. exploit special characteristics of problem
+ //----------------------------------------------------------------------------
+
+ std::cerr << "exploit detected special properties: ";
+ if(np2->hessian_constant())
+ {
+ std::cerr << "*constant hessian* ";
+ app().Options()->SetStringValue("hessian_constant", "yes");
+ }
+
+ if(np2->jac_c_constant())
+ {
+ std::cerr << "*constant jacobian of equality constraints* ";
+ app().Options()->SetStringValue("jac_c_constant", "yes");
+ }
+
+ if(np2->jac_d_constant())
+ {
+ std::cerr << "*constant jacobian of in-equality constraints*";
+ app().Options()->SetStringValue("jac_d_constant", "yes");
+ }
+ std::cerr << std::endl;
+
+ //----------------------------------------------------------------------------
+ // 3. solve problem
+ //----------------------------------------------------------------------------
+
+ // Initialize the IpoptApplication and process the options
+ Ipopt::ApplicationReturnStatus status;
+ status = app_->Initialize();
+ if (status != Ipopt::Solve_Succeeded)
+ {
+ printf("\n\n*** Error IPOPT during initialization!\n");
+ }
+
+ status = app_->OptimizeTNLP( np);
+
+ //----------------------------------------------------------------------------
+ // 4. output statistics
+ //----------------------------------------------------------------------------
+ if (status == Ipopt::Solve_Succeeded || status == Ipopt::Solved_To_Acceptable_Level)
+ {
+ // Retrieve some statistics about the solve
+ Ipopt::Index iter_count = app_->Statistics()->IterationCount();
+ printf("\n\n*** IPOPT: The problem solved in %d iterations!\n", iter_count);
+
+ Ipopt::Number final_obj = app_->Statistics()->FinalObjective();
+ printf("\n\n*** IPOPT: The final value of the objective function is %e.\n", final_obj);
+ }
+
+ return status;
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+
+int
+IPOPTSolver::
+solve(NProblemInterface* _problem,
+ const std::vector& _constraints,
+ const std::vector& _lazy_constraints,
+ const double _almost_infeasible,
+ const int _max_passes )
+{
+ //----------------------------------------------------------------------------
+ // 0. Initialize IPOPT Applicaiton
+ //----------------------------------------------------------------------------
+
+ StopWatch sw; sw.start();
+
+ // Initialize the IpoptApplication and process the options
+ Ipopt::ApplicationReturnStatus status;
+ status = app_->Initialize();
+ if (status != Ipopt::Solve_Succeeded)
+ {
+ printf("\n\n*** Error IPOPT during initialization!\n");
+ }
+
+ bool feasible_point_found = false;
+ int cur_pass = 0;
+ double acceptable_tolerance = 0.01; // hack: read out from ipopt!!!
+ // copy default constraints
+ std::vector constraints = _constraints;
+ std::vector lazy_added(_lazy_constraints.size(),false);
+
+ // cache statistics of all iterations
+ std::vector n_inf;
+ std::vector n_almost_inf;
+
+ while(!feasible_point_found && cur_pass <(_max_passes-1))
+ {
+ ++cur_pass;
+ //----------------------------------------------------------------------------
+ // 1. Create an instance of current IPOPT NLP
+ //----------------------------------------------------------------------------
+ Ipopt::SmartPtr np = new NProblemIPOPT(_problem, constraints);
+ NProblemIPOPT* np2 = dynamic_cast (Ipopt::GetRawPtr(np));
+ // enable caching of solution
+ np2->store_solution() = true;
+
+ //----------------------------------------------------------------------------
+ // 2. exploit special characteristics of problem
+ //----------------------------------------------------------------------------
+
+ std::cerr << "detected special properties which will be exploit: ";
+ if(np2->hessian_constant())
+ {
+ std::cerr << "*constant hessian* ";
+ app().Options()->SetStringValue("hessian_constant", "yes");
+ }
+
+ if(np2->jac_c_constant())
+ {
+ std::cerr << "*constant jacobian of equality constraints* ";
+ app().Options()->SetStringValue("jac_c_constant", "yes");
+ }
+
+ if(np2->jac_d_constant())
+ {
+ std::cerr << "*constant jacobian of in-equality constraints*";
+ app().Options()->SetStringValue("jac_d_constant", "yes");
+ }
+ std::cerr << std::endl;
+
+ //----------------------------------------------------------------------------
+ // 3. solve problem
+ //----------------------------------------------------------------------------
+ status = app_->OptimizeTNLP( np);
+
+ // check lazy constraints
+ n_inf.push_back(0);
+ n_almost_inf.push_back(0);
+ feasible_point_found = true;
+ for(unsigned int i=0; i<_lazy_constraints.size(); ++i)
+ if(!lazy_added[i])
+ {
+ NConstraintInterface* lc = _lazy_constraints[i];
+
+ double v = lc->eval_constraint(&(np2->solution()[0]));
+
+ bool inf = false;
+ bool almost_inf = false;
+
+ if(lc->constraint_type() == NConstraintInterface::NC_EQUAL)
+ {
+ v = std::abs(v);
+ if(v>acceptable_tolerance)
+ inf = true;
+ else
+ if(v>_almost_infeasible)
+ almost_inf = true;
+ }
+ else
+ if(lc->constraint_type() == NConstraintInterface::NC_GREATER_EQUAL)
+ {
+ if(v<-acceptable_tolerance)
+ inf = true;
+ else
+ if(v<_almost_infeasible)
+ almost_inf = true;
+ }
+ else
+ if(lc->constraint_type() == NConstraintInterface::NC_LESS_EQUAL)
+ {
+ if(v>acceptable_tolerance)
+ inf = true;
+ else
+ if(v>-_almost_infeasible)
+ almost_inf = true;
+ }
+
+ // infeasible?
+ if(inf)
+ {
+ constraints.push_back(lc);
+ lazy_added[i] = true;
+ feasible_point_found = false;
+ ++n_inf.back();
+ }
+
+ // almost violated or violated? -> add to constraints
+ if(almost_inf)
+ {
+ constraints.push_back(lc);
+ lazy_added[i] = true;
+ ++n_almost_inf.back();
+ }
+ }
+ }
+
+ // no termination after max number of passes?
+ if(!feasible_point_found)
+ {
+ ++cur_pass;
+
+ std::cerr << "*************** could not find feasible point after " << _max_passes-1 << " -> solving with all lazy constraints..." << std::endl;
+ for(unsigned int i=0; i<_lazy_constraints.size(); ++i)
+ if(!lazy_added[i])
+ constraints.push_back(_lazy_constraints[i]);
+
+ //----------------------------------------------------------------------------
+ // 1. Create an instance of current IPOPT NLP
+ //----------------------------------------------------------------------------
+ Ipopt::SmartPtr np = new NProblemIPOPT(_problem, constraints);
+ NProblemIPOPT* np2 = dynamic_cast (Ipopt::GetRawPtr(np));
+ // enable caching of solution
+ np2->store_solution() = true;
+
+ //----------------------------------------------------------------------------
+ // 2. exploit special characteristics of problem
+ //----------------------------------------------------------------------------
+
+ std::cerr << "exploit detected special properties: ";
+ if(np2->hessian_constant())
+ {
+ std::cerr << "*constant hessian* ";
+ app().Options()->SetStringValue("hessian_constant", "yes");
+ }
+
+ if(np2->jac_c_constant())
+ {
+ std::cerr << "*constant jacobian of equality constraints* ";
+ app().Options()->SetStringValue("jac_c_constant", "yes");
+ }
+
+ if(np2->jac_d_constant())
+ {
+ std::cerr << "*constant jacobian of in-equality constraints*";
+ app().Options()->SetStringValue("jac_d_constant", "yes");
+ }
+ std::cerr << std::endl;
+
+ //----------------------------------------------------------------------------
+ // 3. solve problem
+ //----------------------------------------------------------------------------
+ status = app_->OptimizeTNLP( np);
+ }
+
+ const double overall_time = sw.stop()/1000.0;
+
+ //----------------------------------------------------------------------------
+ // 4. output statistics
+ //----------------------------------------------------------------------------
+ if (status == Ipopt::Solve_Succeeded || status == Ipopt::Solved_To_Acceptable_Level)
+ {
+ // Retrieve some statistics about the solve
+ Ipopt::Index iter_count = app_->Statistics()->IterationCount();
+ printf("\n\n*** IPOPT: The problem solved in %d iterations!\n", iter_count);
+
+ Ipopt::Number final_obj = app_->Statistics()->FinalObjective();
+ printf("\n\n*** IPOPT: The final value of the objective function is %e.\n", final_obj);
+ }
+
+ std::cerr <<"############# IPOPT with lazy constraints statistics ###############" << std::endl;
+ std::cerr << "overall time: " << overall_time << "s" << std::endl;
+ std::cerr << "#passes : " << cur_pass << "( of " << _max_passes << ")" << std::endl;
+ for(unsigned int i=0; i constraints;
+ return this->solve(_problem, constraints);
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+int
+IPOPTSolver::
+solve(NProblemGmmInterface* _problem, std::vector& _constraints)
+{
+ std::cerr << "****** Warning: NProblemGmmInterface is deprecated!!! -> use NProblemInterface *******\n";
+
+ //----------------------------------------------------------------------------
+ // 1. Create an instance of IPOPT NLP
+ //----------------------------------------------------------------------------
+ Ipopt::SmartPtr np = new NProblemGmmIPOPT(_problem, _constraints);
+
+ //----------------------------------------------------------------------------
+ // 2. solve problem
+ //----------------------------------------------------------------------------
+
+ // Initialize the IpoptApplication and process the options
+ Ipopt::ApplicationReturnStatus status;
+ status = app_->Initialize();
+ if (status != Ipopt::Solve_Succeeded)
+ {
+ printf("\n\n*** Error IPOPT during initialization!\n");
+ }
+
+ //----------------------------------------------------------------------------
+ // 3. solve problem
+ //----------------------------------------------------------------------------
+ status = app_->OptimizeTNLP(np);
+
+ //----------------------------------------------------------------------------
+ // 4. output statistics
+ //----------------------------------------------------------------------------
+ if (status == Ipopt::Solve_Succeeded || status == Ipopt::Solved_To_Acceptable_Level)
+ {
+ // Retrieve some statistics about the solve
+ Ipopt::Index iter_count = app_->Statistics()->IterationCount();
+ printf("\n\n*** IPOPT: The problem solved in %d iterations!\n", iter_count);
+
+ Ipopt::Number final_obj = app_->Statistics()->FinalObjective();
+ printf("\n\n*** IPOPT: The final value of the objective function is %e.\n", final_obj);
+ }
+
+ return status;
+}
+
+
+//== IMPLEMENTATION PROBLEM INSTANCE==========================================================
+
+
+void
+NProblemIPOPT::
+split_constraints(const std::vector& _constraints)
+{
+ // split user-provided constraints into general-constraints and bound-constraints
+ constraints_ .clear(); constraints_.reserve(_constraints.size());
+ bound_constraints_.clear(); bound_constraints_.reserve(_constraints.size());
+
+ for(unsigned int i=0; i<_constraints.size(); ++i)
+ {
+ BoundConstraint* bnd_ptr = dynamic_cast(_constraints[i]);
+
+ if(bnd_ptr)
+ bound_constraints_.push_back(bnd_ptr);
+ else
+ constraints_.push_back(_constraints[i]);
+ }
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+void
+NProblemIPOPT::
+analyze_special_properties(const NProblemInterface* _problem, const std::vector& _constraints)
+{
+ hessian_constant_ = true;
+ jac_c_constant_ = true;
+ jac_d_constant_ = true;
+
+ if(!_problem->constant_hessian())
+ hessian_constant_ = false;
+
+ for(unsigned int i=0; i<_constraints.size(); ++i)
+ {
+ if(!_constraints[i]->constant_hessian())
+ hessian_constant_ = false;
+
+ if(!_constraints[i]->constant_gradient())
+ {
+ if(_constraints[i]->constraint_type() == NConstraintInterface::NC_EQUAL)
+ jac_c_constant_ = false;
+ else
+ jac_d_constant_ = false;
+ }
+
+ // nothing else to check?
+ if(!hessian_constant_ && !jac_c_constant_ && !jac_d_constant_)
+ break;
+ }
+
+ //hessian of Lagrangian is only constant, if all hessians of the constraints are zero (due to lambda multipliers)
+ if(!jac_c_constant_ || !jac_d_constant_)
+ hessian_constant_ = false;
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+bool NProblemIPOPT::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
+ Index& nnz_h_lag, IndexStyleEnum& index_style)
+{
+ // number of variables
+ n = problem_->n_unknowns();
+
+ // number of constraints
+ m = constraints_.size();
+
+ // get non-zeros of hessian of lagrangian and jacobi of constraints
+ nnz_jac_g = 0;
+ nnz_h_lag = 0;
+
+ // get nonzero structure
+ std::vector x(n);
+ problem_->initial_x(P(x));
+
+ // nonzeros in the jacobian of C_ and the hessian of the lagrangian
+ SMatrixNP HP;
+ SVectorNC g;
+ SMatrixNC H;
+ problem_->eval_hessian(P(x), HP);
+
+ // get nonzero structure of hessian of problem
+ for(int i=0; i= it.col())
+ ++nnz_h_lag;
+
+ // get nonzero structure of constraints
+ for( int i=0; ieval_gradient(P(x),g);
+
+ nnz_jac_g += g.nonZeros();
+
+ // count lower triangular elements
+ constraints_[i]->eval_hessian (P(x),H);
+
+ SMatrixNC::iterator m_it = H.begin();
+ for(; m_it != H.end(); ++m_it)
+ if( m_it.row() >= m_it.col())
+ ++nnz_h_lag;
+ }
+
+ // We use the standard fortran index style for row/col entries
+ index_style = C_STYLE;
+
+ return true;
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+bool NProblemIPOPT::get_bounds_info(Index n, Number* x_l, Number* x_u,
+ Index m, Number* g_l, Number* g_u)
+{
+ // check dimensions
+ if( n != (Index)problem_->n_unknowns())
+ std::cerr << "Warning: IPOPT #unknowns != n " << n << problem_->n_unknowns() << std::endl;
+ if( m != (Index)constraints_.size())
+ std::cerr << "Warning: IPOPT #constraints != m " << m << constraints_.size() << std::endl;
+
+
+ // first clear all variable bounds
+ for( int i=0; iidx()) < n)
+ {
+ switch(bound_constraints_[i]->constraint_type())
+ {
+ case NConstraintInterface::NC_LESS_EQUAL:
+ {
+ x_u[bound_constraints_[i]->idx()] = bound_constraints_[i]->bound();
+ }break;
+
+ case NConstraintInterface::NC_GREATER_EQUAL:
+ {
+ x_l[bound_constraints_[i]->idx()] = bound_constraints_[i]->bound();
+ }break;
+
+ case NConstraintInterface::NC_EQUAL:
+ {
+ x_l[bound_constraints_[i]->idx()] = bound_constraints_[i]->bound();
+ x_u[bound_constraints_[i]->idx()] = bound_constraints_[i]->bound();
+ }break;
+ }
+ }
+ else
+ std::cerr << "Warning: invalid bound constraint in IPOPTSolver!!!" << std::endl;
+ }
+
+ // set bounds for constraints
+ for( int i=0; iconstraint_type())
+ {
+ case NConstraintInterface::NC_EQUAL : g_u[i] = 0.0 ; g_l[i] = 0.0 ; break;
+ case NConstraintInterface::NC_LESS_EQUAL : g_u[i] = 0.0 ; g_l[i] = -1.0e19; break;
+ case NConstraintInterface::NC_GREATER_EQUAL : g_u[i] = 1.0e19; g_l[i] = 0.0 ; break;
+ default : g_u[i] = 1.0e19; g_l[i] = -1.0e19; break;
+ }
+ }
+
+ return true;
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+bool NProblemIPOPT::get_starting_point(Index n, bool init_x, Number* x,
+ bool init_z, Number* z_L, Number* z_U,
+ Index m, bool init_lambda,
+ Number* lambda)
+{
+ // get initial value of problem instance
+ problem_->initial_x(x);
+
+ return true;
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+bool NProblemIPOPT::eval_f(Index n, const Number* x, bool new_x, Number& obj_value)
+{
+ // return the value of the objective function
+ obj_value = problem_->eval_f(x);
+ return true;
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+bool NProblemIPOPT::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f)
+{
+ problem_->eval_gradient(x, grad_f);
+
+ return true;
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+bool NProblemIPOPT::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
+{
+ // evaluate all constraint functions
+ for( int i=0; ieval_constraint(x);
+
+ return true;
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+bool NProblemIPOPT::eval_jac_g(Index n, const Number* x, bool new_x,
+ Index m, Index nele_jac, Index* iRow, Index *jCol,
+ Number* values)
+{
+ if (values == NULL)
+ {
+ // get x for evaluation (arbitrary position should be ok)
+ std::vector x_rnd(problem_->n_unknowns(), 0.0);
+
+ int gi = 0;
+ SVectorNC g;
+ for( int i=0; ieval_gradient(&(x_rnd[0]), g);
+ SVectorNC::InnerIterator v_it(g);
+ for( ; v_it; ++v_it)
+ {
+ iRow[gi] = i;
+ jCol[gi] = v_it.index();
+ ++gi;
+ }
+ }
+ }
+ else
+ {
+ // return the values of the jacobian of the constraints
+
+ // return the structure of the jacobian of the constraints
+ // global index
+ int gi = 0;
+ SVectorNC g;
+
+ for( int i=0; ieval_gradient(x, g);
+
+ SVectorNC::InnerIterator v_it(g);
+
+ for( ; v_it; ++v_it)
+ {
+ values[gi] = v_it.value();
+ ++gi;
+ }
+ }
+
+ if( gi != nele_jac)
+ std::cerr << "Warning: number of non-zeros in Jacobian of C is incorrect: "
+ << gi << " vs " << nele_jac << std::endl;
+ }
+
+ return true;
+}
+
+
+//-----------------------------------------------------------------------------
+
+
+bool NProblemIPOPT::eval_h(Index n, const Number* x, bool new_x,
+ Number obj_factor, Index m, const Number* lambda,
+ bool new_lambda, Index nele_hess, Index* iRow,
+ Index* jCol, Number* values)
+{
+ if (values == NULL)
+ {
+ // return structure
+
+ // get x for evaluation (arbitrary position should be ok)
+ std::vector x_rnd(problem_->n_unknowns(), 0.0);
+
+ // global index
+ int gi = 0;
+ // get hessian of problem
+ SMatrixNP HP;
+ problem_->eval_hessian(&(x_rnd[0]), HP);
+
+ for(int i=0; i= it.col())
+ {
+ // it.value();
+ iRow[gi] = it.row();
+ jCol[gi] = it.col();
+ ++gi;
+ }
+ }
+
+ // Hessians of Constraints
+ for(unsigned int j=0; jeval_hessian(&(x_rnd[0]), H);
+
+ SMatrixNC::iterator m_it = H.begin();
+ SMatrixNC::iterator m_end = H.end();
+
+ for(; m_it != m_end; ++m_it)
+ {
+ // store lower triangular part only
+ if( m_it.row() >= m_it.col())
+ {
+ iRow[gi] = m_it.row();
+ jCol[gi] = m_it.col();
+ ++gi;
+ }
+ }
+ }
+
+ // error check
+ if( gi != nele_hess)
+ std::cerr << "Warning: number of non-zeros in Hessian of Lagrangian is incorrect while indexing: "
+ << gi << " vs " << nele_hess << std::endl;
+ }
+ else
+ {
+ // return values.
+
+ // global index
+ int gi = 0;
+ // get hessian of problem
+ SMatrixNP HP;
+ problem_->eval_hessian(x, HP);
+
+ for(int i=0; i