diff options
Diffstat (limited to 'xs/src/igl/mosek')
| -rw-r--r-- | xs/src/igl/mosek/bbw.cpp | 88 | ||||
| -rw-r--r-- | xs/src/igl/mosek/bbw.h | 60 | ||||
| -rw-r--r-- | xs/src/igl/mosek/mosek_guarded.cpp | 24 | ||||
| -rw-r--r-- | xs/src/igl/mosek/mosek_guarded.h | 31 | ||||
| -rw-r--r-- | xs/src/igl/mosek/mosek_linprog.cpp | 164 | ||||
| -rw-r--r-- | xs/src/igl/mosek/mosek_linprog.h | 59 | ||||
| -rw-r--r-- | xs/src/igl/mosek/mosek_quadprog.cpp | 343 | ||||
| -rw-r--r-- | xs/src/igl/mosek/mosek_quadprog.h | 145 |
8 files changed, 0 insertions, 914 deletions
diff --git a/xs/src/igl/mosek/bbw.cpp b/xs/src/igl/mosek/bbw.cpp deleted file mode 100644 index 99d28ab15..000000000 --- a/xs/src/igl/mosek/bbw.cpp +++ /dev/null @@ -1,88 +0,0 @@ -// This file is part of libigl, a simple c++ geometry processing library. -// -// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com> -// -// This Source Code Form is subject to the terms of the Mozilla Public License -// v. 2.0. If a copy of the MPL was not distributed with this file, You can -// obtain one at http://mozilla.org/MPL/2.0/. -#include "bbw.h" -#include "mosek_quadprog.h" -#include "../harmonic.h" -#include "../slice_into.h" -#include <Eigen/Sparse> -#include <iostream> -#include <cstdio> - - -template < - typename DerivedV, - typename DerivedEle, - typename Derivedb, - typename Derivedbc, - typename DerivedW> -IGL_INLINE bool igl::mosek::bbw( - const Eigen::PlainObjectBase<DerivedV> & V, - const Eigen::PlainObjectBase<DerivedEle> & Ele, - const Eigen::PlainObjectBase<Derivedb> & b, - const Eigen::PlainObjectBase<Derivedbc> & bc, - igl::BBWData & data, - igl::mosek::MosekData & mosek_data, - Eigen::PlainObjectBase<DerivedW> & W - ) -{ - using namespace std; - using namespace Eigen; - assert(!data.partition_unity && "partition_unity not implemented yet"); - // number of domain vertices - int n = V.rows(); - // number of handles - int m = bc.cols(); - // Build biharmonic operator - Eigen::SparseMatrix<typename DerivedV::Scalar> Q; - harmonic(V,Ele,2,Q); - W.derived().resize(n,m); - // No linear terms - VectorXd c = VectorXd::Zero(n); - // No linear constraints - SparseMatrix<typename DerivedW::Scalar> A(0,n); - VectorXd uc(0,1),lc(0,1); - // Upper and lower box constraints (Constant bounds) - VectorXd ux = VectorXd::Ones(n); - VectorXd lx = VectorXd::Zero(n); - // Loop over handles - for(int i = 0;i<m;i++) - { - if(data.verbosity >= 1) - { - cout<<"BBW: Computing weight for handle "<<i+1<<" out of "<<m<< - "."<<endl; - } - VectorXd bci = bc.col(i); - VectorXd Wi; - // impose boundary conditions via bounds - slice_into(bci,b,ux); - slice_into(bci,b,lx); - bool r = mosek_quadprog(Q,c,0,A,lc,uc,lx,ux,mosek_data,Wi); - if(!r) - { - return false; - } - W.col(i) = Wi; - } -#ifndef NDEBUG - const double min_rowsum = W.rowwise().sum().array().abs().minCoeff(); - if(min_rowsum < 0.1) - { - cerr<<"bbw.cpp: Warning, minimum row sum is very low. Consider more " - "active set iterations or enforcing partition of unity."<<endl; - } -#endif - - return true; -} - -#ifdef IGL_STATIC_LIBRARY -// Explicit template instantiation -template bool igl::mosek::bbw<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, igl::BBWData&, igl::mosek::MosekData&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&); -#endif - diff --git a/xs/src/igl/mosek/bbw.h b/xs/src/igl/mosek/bbw.h deleted file mode 100644 index 5dcae0fed..000000000 --- a/xs/src/igl/mosek/bbw.h +++ /dev/null @@ -1,60 +0,0 @@ -// This file is part of libigl, a simple c++ geometry processing library. -// -// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com> -// -// This Source Code Form is subject to the terms of the Mozilla Public License -// v. 2.0. If a copy of the MPL was not distributed with this file, You can -// obtain one at http://mozilla.org/MPL/2.0/. -#ifndef IGL_MOSEK_BBW_H -#define IGL_MOSEK_BBW_H -#include "../igl_inline.h" -#include "mosek_quadprog.h" -#include "../bbw.h" -#include <Eigen/Dense> - -namespace igl -{ - namespace mosek - { - // Compute Bounded Biharmonic Weights on a given domain (V,Ele) with a given - // set of boundary conditions - // - // Templates - // DerivedV derived type of eigen matrix for V (e.g. MatrixXd) - // DerivedF derived type of eigen matrix for F (e.g. MatrixXi) - // Derivedb derived type of eigen matrix for b (e.g. VectorXi) - // Derivedbc derived type of eigen matrix for bc (e.g. MatrixXd) - // DerivedW derived type of eigen matrix for W (e.g. MatrixXd) - // Inputs: - // V #V by dim vertex positions - // Ele #Elements by simplex-size list of element indices - // b #b boundary indices into V - // bc #b by #W list of boundary values - // data object containing options, initial guess --> solution and results - // mosek_data object containing mosek options - // Outputs: - // W #V by #W list of *unnormalized* weights to normalize use - // igl::normalize_row_sums(W,W); - // Returns true on success, false on failure - template < - typename DerivedV, - typename DerivedEle, - typename Derivedb, - typename Derivedbc, - typename DerivedW> - IGL_INLINE bool bbw( - const Eigen::PlainObjectBase<DerivedV> & V, - const Eigen::PlainObjectBase<DerivedEle> & Ele, - const Eigen::PlainObjectBase<Derivedb> & b, - const Eigen::PlainObjectBase<Derivedbc> & bc, - igl::BBWData & data, - igl::mosek::MosekData & mosek_data, - Eigen::PlainObjectBase<DerivedW> & W); - } -} - -#ifndef IGL_STATIC_LIBRARY -# include "bbw.cpp" -#endif - -#endif diff --git a/xs/src/igl/mosek/mosek_guarded.cpp b/xs/src/igl/mosek/mosek_guarded.cpp deleted file mode 100644 index 242d928f6..000000000 --- a/xs/src/igl/mosek/mosek_guarded.cpp +++ /dev/null @@ -1,24 +0,0 @@ -// This file is part of libigl, a simple c++ geometry processing library. -// -// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com> -// -// This Source Code Form is subject to the terms of the Mozilla Public License -// v. 2.0. If a copy of the MPL was not distributed with this file, You can -// obtain one at http://mozilla.org/MPL/2.0/. -#include "mosek_guarded.h" -#include <iostream> - -IGL_INLINE MSKrescodee igl::mosek::mosek_guarded(const MSKrescodee r) -{ - using namespace std; - if(r != MSK_RES_OK) - { - /* In case of an error print error code and description. */ - char symname[MSK_MAX_STR_LEN]; - char desc[MSK_MAX_STR_LEN]; - MSK_getcodedesc(r,symname,desc); - cerr<<"MOSEK ERROR ("<<r<<"): "<<symname<<" - '"<<desc<<"'"<<endl; - } - return r; -} - diff --git a/xs/src/igl/mosek/mosek_guarded.h b/xs/src/igl/mosek/mosek_guarded.h deleted file mode 100644 index 1c9e0cbf9..000000000 --- a/xs/src/igl/mosek/mosek_guarded.h +++ /dev/null @@ -1,31 +0,0 @@ -// This file is part of libigl, a simple c++ geometry processing library. -// -// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com> -// -// This Source Code Form is subject to the terms of the Mozilla Public License -// v. 2.0. If a copy of the MPL was not distributed with this file, You can -// obtain one at http://mozilla.org/MPL/2.0/. -#ifndef IGL_MOSEK_MOSEK_GUARDED_H -#define IGL_MOSEK_MOSEK_GUARDED_H -#include "../igl_inline.h" - -#include "mosek.h" -namespace igl -{ - namespace mosek - { - // Little function to wrap around mosek call to handle errors - // - // Inputs: - // r mosek error code returned from mosek call - // Returns r untouched - IGL_INLINE MSKrescodee mosek_guarded(const MSKrescodee r); - } -} - -#ifndef IGL_STATIC_LIBRARY -# include "mosek_guarded.cpp" -#endif - -#endif - diff --git a/xs/src/igl/mosek/mosek_linprog.cpp b/xs/src/igl/mosek/mosek_linprog.cpp deleted file mode 100644 index fe09d5adb..000000000 --- a/xs/src/igl/mosek/mosek_linprog.cpp +++ /dev/null @@ -1,164 +0,0 @@ -// This file is part of libigl, a simple c++ geometry processing library. -// -// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com> -// -// This Source Code Form is subject to the terms of the Mozilla Public License -// v. 2.0. If a copy of the MPL was not distributed with this file, You can -// obtain one at http://mozilla.org/MPL/2.0/. -#include "mosek_linprog.h" -#include "../mosek/mosek_guarded.h" -#include "../harwell_boeing.h" -#include <limits> -#include <cmath> -#include <vector> - -IGL_INLINE bool igl::mosek::mosek_linprog( - const Eigen::VectorXd & c, - const Eigen::SparseMatrix<double> & A, - const Eigen::VectorXd & lc, - const Eigen::VectorXd & uc, - const Eigen::VectorXd & lx, - const Eigen::VectorXd & ux, - Eigen::VectorXd & x) -{ - // variables for mosek task, env and result code - MSKenv_t env; - // Create the MOSEK environment - mosek_guarded(MSK_makeenv(&env,NULL)); - // initialize mosek environment -#if MSK_VERSION_MAJOR <= 7 - mosek_guarded(MSK_initenv(env)); -#endif - const bool ret = mosek_linprog(c,A,lc,uc,lx,ux,env,x); - MSK_deleteenv(&env); - return ret; -} - -IGL_INLINE bool igl::mosek::mosek_linprog( - const Eigen::VectorXd & c, - const Eigen::SparseMatrix<double> & A, - const Eigen::VectorXd & lc, - const Eigen::VectorXd & uc, - const Eigen::VectorXd & lx, - const Eigen::VectorXd & ux, - const MSKenv_t & env, - Eigen::VectorXd & x) -{ - // following http://docs.mosek.com/7.1/capi/Linear_optimization.html - using namespace std; - // number of constraints - const int m = A.rows(); - // number of variables - const int n = A.cols(); - - - vector<double> vAv; - vector<int> vAri,vAcp; - int nr; - harwell_boeing(A,nr,vAv,vAri,vAcp); - - MSKtask_t task; - // Create the optimization task - mosek_guarded(MSK_maketask(env,m,n,&task)); - // no threads - mosek_guarded(MSK_putintparam(task,MSK_IPAR_NUM_THREADS,1)); - if(m>0) - { - // Append 'm' empty constraints, the constrainst will initially have no - // bounds - mosek_guarded(MSK_appendcons(task,m)); - } - mosek_guarded(MSK_appendvars(task,n)); - - - const auto & key = [](const double lxj, const double uxj) -> - MSKboundkeye - { - MSKboundkeye k = MSK_BK_FR; - if(isfinite(lxj) && isfinite(uxj)) - { - if(lxj == uxj) - { - k = MSK_BK_FX; - }else{ - k = MSK_BK_RA; - } - }else if(isfinite(lxj)) - { - k = MSK_BK_LO; - }else if(isfinite(uxj)) - { - k = MSK_BK_UP; - } - return k; - }; - - // loop over variables - for(int j = 0;j<n;j++) - { - if(c.size() > 0) - { - // Set linear term c_j in the objective - mosek_guarded(MSK_putcj(task,j,c(j))); - } - - // Set constant bounds on variable j - const double lxj = lx.size()>0?lx[j]:-numeric_limits<double>::infinity(); - const double uxj = ux.size()>0?ux[j]: numeric_limits<double>::infinity(); - mosek_guarded(MSK_putvarbound(task,j,key(lxj,uxj),lxj,uxj)); - - if(m>0) - { - // Input column j of A - mosek_guarded( - MSK_putacol( - task, - j, - vAcp[j+1]-vAcp[j], - &vAri[vAcp[j]], - &vAv[vAcp[j]]) - ); - } - } - // loop over constraints - for(int i = 0;i<m;i++) - { - // Set constraint bounds for row i - const double lci = lc.size()>0?lc[i]:-numeric_limits<double>::infinity(); - const double uci = uc.size()>0?uc[i]: numeric_limits<double>::infinity(); - mosek_guarded(MSK_putconbound(task,i,key(lci,uci),lci,uci)); - } - - // Now the optimizer has been prepared - MSKrescodee trmcode; - // run the optimizer - mosek_guarded(MSK_optimizetrm(task,&trmcode)); - // Get status - MSKsolstae solsta; - MSK_getsolsta (task,MSK_SOL_ITR,&solsta); - bool success = false; - switch(solsta) - { - case MSK_SOL_STA_OPTIMAL: - case MSK_SOL_STA_NEAR_OPTIMAL: - x.resize(n); - /* Request the basic solution. */ - MSK_getxx(task,MSK_SOL_BAS,x.data()); - success = true; - break; - case MSK_SOL_STA_DUAL_INFEAS_CER: - case MSK_SOL_STA_PRIM_INFEAS_CER: - case MSK_SOL_STA_NEAR_DUAL_INFEAS_CER: - case MSK_SOL_STA_NEAR_PRIM_INFEAS_CER: - //printf("Primal or dual infeasibility certificate found.\n"); - break; - case MSK_SOL_STA_UNKNOWN: - //printf("The status of the solution could not be determined.\n"); - break; - default: - //printf("Other solution status."); - break; - } - MSK_deletetask(&task); - return success; -} diff --git a/xs/src/igl/mosek/mosek_linprog.h b/xs/src/igl/mosek/mosek_linprog.h deleted file mode 100644 index 57791cd60..000000000 --- a/xs/src/igl/mosek/mosek_linprog.h +++ /dev/null @@ -1,59 +0,0 @@ -// This file is part of libigl, a simple c++ geometry processing library. -// -// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com> -// -// This Source Code Form is subject to the terms of the Mozilla Public License -// v. 2.0. If a copy of the MPL was not distributed with this file, You can -// obtain one at http://mozilla.org/MPL/2.0/. -#ifndef IGL_MOSEK_MOSEK_LINPROG_H -#define IGL_MOSEK_MOSEK_LINPROG_H -#include "../igl_inline.h" -#include <Eigen/Core> -#include <Eigen/Sparse> -#include <mosek.h> -namespace igl -{ - namespace mosek - { - // Solve a linear program using mosek: - // - // min c'x - // s.t. lc <= A x <= uc - // lx <= x <= ux - // - // Inputs: - // c #x list of linear objective coefficients - // A #A by #x matrix of linear inequality constraint coefficients - // lc #A list of lower constraint bounds - // uc #A list of upper constraint bounds - // lx #x list of lower variable bounds - // ux #x list of upper variable bounds - // Outputs: - // x #x list of solution values - // Returns true iff success. - IGL_INLINE bool mosek_linprog( - const Eigen::VectorXd & c, - const Eigen::SparseMatrix<double> & A, - const Eigen::VectorXd & lc, - const Eigen::VectorXd & uc, - const Eigen::VectorXd & lx, - const Eigen::VectorXd & ux, - Eigen::VectorXd & x); - // Wrapper that keeps mosek environment alive (if licence checking is - // becoming a bottleneck) - IGL_INLINE bool mosek_linprog( - const Eigen::VectorXd & c, - const Eigen::SparseMatrix<double> & A, - const Eigen::VectorXd & lc, - const Eigen::VectorXd & uc, - const Eigen::VectorXd & lx, - const Eigen::VectorXd & ux, - const MSKenv_t & env, - Eigen::VectorXd & x); - } -} - -#ifndef IGL_STATIC_LIBRARY -# include "mosek_linprog.cpp" -#endif -#endif diff --git a/xs/src/igl/mosek/mosek_quadprog.cpp b/xs/src/igl/mosek/mosek_quadprog.cpp deleted file mode 100644 index 36d9e873b..000000000 --- a/xs/src/igl/mosek/mosek_quadprog.cpp +++ /dev/null @@ -1,343 +0,0 @@ -// This file is part of libigl, a simple c++ geometry processing library. -// -// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com> -// -// This Source Code Form is subject to the terms of the Mozilla Public License -// v. 2.0. If a copy of the MPL was not distributed with this file, You can -// obtain one at http://mozilla.org/MPL/2.0/. -#include "mosek_quadprog.h" -#include "mosek_guarded.h" -#include <cstdio> -#include "../find.h" -#include "../verbose.h" -#include "../speye.h" -#include "../matrix_to_list.h" -#include "../list_to_matrix.h" -#include "../harwell_boeing.h" -#include "../EPS.h" - - -igl::mosek::MosekData::MosekData() -{ - // These are the default settings that worked well for BBW. Your miles may - // very well be kilometers. - - // >1e0 NONSOLUTION - // 1e-1 artifacts in deformation - // 1e-3 artifacts in isolines - // 1e-4 seems safe - // 1e-8 MOSEK DEFAULT SOLUTION - douparam[MSK_DPAR_INTPNT_TOL_REL_GAP]=1e-8; -#if MSK_VERSION_MAJOR >= 8 - douparam[MSK_DPAR_INTPNT_QO_TOL_REL_GAP]=1e-12; -#endif - // Force using multiple threads, not sure if MOSEK is properly destroying - //extra threads... -#if MSK_VERSION_MAJOR >= 7 - intparam[MSK_IPAR_NUM_THREADS] = 6; -#elif MSK_VERSION_MAJOR == 6 - intparam[MSK_IPAR_INTPNT_NUM_THREADS] = 6; -#endif -#if MSK_VERSION_MAJOR == 6 - // Force turn off data check - intparam[MSK_IPAR_DATA_CHECK]=MSK_OFF; -#endif - // Turn off presolving - // intparam[MSK_IPAR_PRESOLVE_USE] = MSK_PRESOLVE_MODE_OFF; - // Force particular matrix reordering method - // MSK_ORDER_METHOD_NONE cuts time in half roughly, since half the time is - // usually spent reordering the matrix - // !! WARNING Setting this parameter to anything but MSK_ORDER_METHOD_FREE - // seems to have the effect of setting it to MSK_ORDER_METHOD_NONE - // *Or maybe Mosek is spending a bunch of time analyzing the matrix to - // choose the right ordering method when really any of them are - // instantaneous - intparam[MSK_IPAR_INTPNT_ORDER_METHOD] = MSK_ORDER_METHOD_NONE; - // 1.0 means optimizer is very lenient about declaring model infeasible - douparam[MSK_DPAR_INTPNT_TOL_INFEAS] = 1e-8; - // Hard to say if this is doing anything, probably nothing dramatic - douparam[MSK_DPAR_INTPNT_TOL_PSAFE]= 1e2; - // Turn off convexity check - intparam[MSK_IPAR_CHECK_CONVEXITY] = MSK_CHECK_CONVEXITY_NONE; -} - -template <typename Index, typename Scalar> -IGL_INLINE bool igl::mosek::mosek_quadprog( - const Index n, - std::vector<Index> & Qi, - std::vector<Index> & Qj, - std::vector<Scalar> & Qv, - const std::vector<Scalar> & c, - const Scalar cf, - const Index m, - std::vector<Scalar> & Av, - std::vector<Index> & Ari, - const std::vector<Index> & Acp, - const std::vector<Scalar> & lc, - const std::vector<Scalar> & uc, - const std::vector<Scalar> & lx, - const std::vector<Scalar> & ux, - MosekData & mosek_data, - std::vector<Scalar> & x) -{ - // I J V vectors of Q should all be same length - assert(Qv.size() == Qi.size()); - assert(Qv.size() == Qj.size()); - // number of columns in linear constraint matrix must be ≤ number of - // variables - assert( (int)Acp.size() == (n+1)); - // linear bound vectors must be size of number of constraints or empty - assert( ((int)lc.size() == m) || ((int)lc.size() == 0)); - assert( ((int)uc.size() == m) || ((int)uc.size() == 0)); - // constant bound vectors must be size of number of variables or empty - assert( ((int)lx.size() == n) || ((int)lx.size() == 0)); - assert( ((int)ux.size() == n) || ((int)ux.size() == 0)); - - // allocate space for solution in x - x.resize(n); - - // variables for mosek task, env and result code - MSKenv_t env; - MSKtask_t task; - - // Create the MOSEK environment -#if MSK_VERSION_MAJOR >= 7 - mosek_guarded(MSK_makeenv(&env,NULL)); -#elif MSK_VERSION_MAJOR == 6 - mosek_guarded(MSK_makeenv(&env,NULL,NULL,NULL,NULL)); -#endif - ///* Directs the log stream to the 'printstr' function. */ - //// Little function mosek needs in order to know how to print to std out - //const auto & printstr = [](void *handle, char str[]) - //{ - // printf("%s",str); - //} - //mosek_guarded(MSK_linkfunctoenvstream(env,MSK_STREAM_LOG,NULL,printstr)); - // initialize mosek environment -#if MSK_VERSION_MAJOR <= 7 - mosek_guarded(MSK_initenv(env)); -#endif - // Create the optimization task - mosek_guarded(MSK_maketask(env,m,n,&task)); - verbose("Creating task with %ld linear constraints and %ld variables...\n",m,n); - //// Tell mosek how to print to std out - //mosek_guarded(MSK_linkfunctotaskstream(task,MSK_STREAM_LOG,NULL,printstr)); - // Give estimate of number of variables - mosek_guarded(MSK_putmaxnumvar(task,n)); - if(m>0) - { - // Give estimate of number of constraints - mosek_guarded(MSK_putmaxnumcon(task,m)); - // Give estimate of number of non zeros in A - mosek_guarded(MSK_putmaxnumanz(task,Av.size())); - } - // Give estimate of number of non zeros in Q - mosek_guarded(MSK_putmaxnumqnz(task,Qv.size())); - if(m>0) - { - // Append 'm' empty constraints, the constrainst will initially have no - // bounds -#if MSK_VERSION_MAJOR >= 7 - mosek_guarded(MSK_appendcons(task,m)); -#elif MSK_VERSION_MAJOR == 6 - mosek_guarded(MSK_append(task,MSK_ACC_CON,m)); -#endif - } - // Append 'n' variables -#if MSK_VERSION_MAJOR >= 7 - mosek_guarded(MSK_appendvars(task,n)); -#elif MSK_VERSION_MAJOR == 6 - mosek_guarded(MSK_append(task,MSK_ACC_VAR,n)); -#endif - // add a contant term to the objective - mosek_guarded(MSK_putcfix(task,cf)); - - // loop over variables - for(int j = 0;j<n;j++) - { - if(c.size() > 0) - { - // Set linear term c_j in the objective - mosek_guarded(MSK_putcj(task,j,c[j])); - } - - // Set constant bounds on variable j - if(lx[j] == ux[j]) - { - mosek_guarded(MSK_putbound(task,MSK_ACC_VAR,j,MSK_BK_FX,lx[j],ux[j])); - }else - { - mosek_guarded(MSK_putbound(task,MSK_ACC_VAR,j,MSK_BK_RA,lx[j],ux[j])); - } - - if(m>0) - { - // Input column j of A -#if MSK_VERSION_MAJOR >= 7 - mosek_guarded( - MSK_putacol( - task, - j, - Acp[j+1]-Acp[j], - &Ari[Acp[j]], - &Av[Acp[j]]) - ); -#elif MSK_VERSION_MAJOR == 6 - mosek_guarded( - MSK_putavec( - task, - MSK_ACC_VAR, - j, - Acp[j+1]-Acp[j], - &Ari[Acp[j]], - &Av[Acp[j]]) - ); -#endif - } - } - - // loop over constraints - for(int i = 0;i<m;i++) - { - // put bounds on constraints - mosek_guarded(MSK_putbound(task,MSK_ACC_CON,i,MSK_BK_RA,lc[i],uc[i])); - } - - // Input Q for the objective (REMEMBER Q SHOULD ONLY BE LOWER TRIANGLE) - mosek_guarded(MSK_putqobj(task,Qv.size(),&Qi[0],&Qj[0],&Qv[0])); - - // Set up task parameters - for( - std::map<MSKiparame,int>::iterator pit = mosek_data.intparam.begin(); - pit != mosek_data.intparam.end(); - pit++) - { - mosek_guarded(MSK_putintparam(task,pit->first,pit->second)); - } - for( - std::map<MSKdparame,double>::iterator pit = mosek_data.douparam.begin(); - pit != mosek_data.douparam.end(); - pit++) - { - mosek_guarded(MSK_putdouparam(task,pit->first,pit->second)); - } - - // Now the optimizer has been prepared - MSKrescodee trmcode; - // run the optimizer - mosek_guarded(MSK_optimizetrm(task,&trmcode)); - - //// Print a summary containing information about the solution for debugging - //// purposes - //MSK_solutionsummary(task,MSK_STREAM_LOG); - - // Get status of solution - MSKsolstae solsta; -#if MSK_VERSION_MAJOR >= 7 - MSK_getsolsta (task,MSK_SOL_ITR,&solsta); -#elif MSK_VERSION_MAJOR == 6 - MSK_getsolutionstatus(task,MSK_SOL_ITR,NULL,&solsta); -#endif - - bool success = false; - switch(solsta) - { - case MSK_SOL_STA_OPTIMAL: - case MSK_SOL_STA_NEAR_OPTIMAL: - MSK_getsolutionslice(task,MSK_SOL_ITR,MSK_SOL_ITEM_XX,0,n,&x[0]); - //printf("Optimal primal solution\n"); - //for(size_t j=0; j<n; ++j) - //{ - // printf("x[%ld]: %g\n",j,x[j]); - //} - success = true; - break; - case MSK_SOL_STA_DUAL_INFEAS_CER: - case MSK_SOL_STA_PRIM_INFEAS_CER: - case MSK_SOL_STA_NEAR_DUAL_INFEAS_CER: - case MSK_SOL_STA_NEAR_PRIM_INFEAS_CER: - //printf("Primal or dual infeasibility certificate found.\n"); - break; - case MSK_SOL_STA_UNKNOWN: - //printf("The status of the solution could not be determined.\n"); - break; - default: - //printf("Other solution status."); - break; - } - - MSK_deletetask(&task); - MSK_deleteenv(&env); - - return success; -} - -IGL_INLINE bool igl::mosek::mosek_quadprog( - const Eigen::SparseMatrix<double> & Q, - const Eigen::VectorXd & c, - const double cf, - const Eigen::SparseMatrix<double> & A, - const Eigen::VectorXd & lc, - const Eigen::VectorXd & uc, - const Eigen::VectorXd & lx, - const Eigen::VectorXd & ux, - MosekData & mosek_data, - Eigen::VectorXd & x) -{ - using namespace Eigen; - using namespace std; - - typedef int Index; - typedef double Scalar; - // Q should be square - assert(Q.rows() == Q.cols()); - // Q should be symmetric -#ifdef EIGEN_HAS_A_BUG_AND_FAILS_TO_LET_ME_COMPUTE_Q_MINUS_Q_TRANSPOSE - assert( (Q-Q.transpose()).sum() < FLOAT_EPS); -#endif - // Only keep lower triangular part of Q - SparseMatrix<Scalar> QL; - //QL = Q.template triangularView<Lower>(); - QL = Q.triangularView<Lower>(); - VectorXi Qi,Qj; - VectorXd Qv; - find(QL,Qi,Qj,Qv); - vector<Index> vQi = matrix_to_list(Qi); - vector<Index> vQj = matrix_to_list(Qj); - vector<Scalar> vQv = matrix_to_list(Qv); - - // Convert linear term - vector<Scalar> vc = matrix_to_list(c); - - assert(lc.size() == A.rows()); - assert(uc.size() == A.rows()); - // Convert A to harwell boeing format - vector<Scalar> vAv; - vector<Index> vAr,vAc; - Index nr; - harwell_boeing(A,nr,vAv,vAr,vAc); - - assert(lx.size() == Q.rows()); - assert(ux.size() == Q.rows()); - vector<Scalar> vlc = matrix_to_list(lc); - vector<Scalar> vuc = matrix_to_list(uc); - vector<Scalar> vlx = matrix_to_list(lx); - vector<Scalar> vux = matrix_to_list(ux); - - vector<Scalar> vx; - bool ret = mosek_quadprog<Index,Scalar>( - Q.rows(),vQi,vQj,vQv, - vc, - cf, - nr, - vAv, vAr, vAc, - vlc,vuc, - vlx,vux, - mosek_data, - vx); - list_to_matrix(vx,x); - return ret; -} -#ifdef IGL_STATIC_LIBRARY -// Explicit template declarations -#endif diff --git a/xs/src/igl/mosek/mosek_quadprog.h b/xs/src/igl/mosek/mosek_quadprog.h deleted file mode 100644 index 38343318a..000000000 --- a/xs/src/igl/mosek/mosek_quadprog.h +++ /dev/null @@ -1,145 +0,0 @@ -// This file is part of libigl, a simple c++ geometry processing library. -// -// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com> -// -// This Source Code Form is subject to the terms of the Mozilla Public License -// v. 2.0. If a copy of the MPL was not distributed with this file, You can -// obtain one at http://mozilla.org/MPL/2.0/. -#ifndef IGL_MOSEK_MOSEK_QUADPROG_H -#define IGL_MOSEK_MOSEK_QUADPROG_H -#include "../igl_inline.h" -#include <vector> -#include <map> -#include <mosek.h> - - -#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET -#include <Eigen/Dense> -#include <Eigen/Sparse> - -namespace igl -{ - namespace mosek - { - struct MosekData - { - // Integer parameters - std::map<MSKiparame,int> intparam; - // Double parameters - std::map<MSKdparame,double> douparam; - // Default values - IGL_INLINE MosekData(); - }; - // Solve a convex quadratic optimization problem with linear and constant - // bounds, that is: - // - // Minimize: ½ * xT * Q⁰ * x + cT * x + cf - // - // Subject to: lc ≤ Ax ≤ uc - // lx ≤ x ≤ ux - // - // where we are trying to find the optimal vector of values x. - // - // Note: Q⁰ must be symmetric and the ½ is a convention of MOSEK - // - // Note: Because of how MOSEK accepts different parts of the system, Q - // should be stored in IJV (aka Coordinate) format and should only include - // entries in the lower triangle. A should be stored in Column compressed - // (aka Harwell Boeing) format. As described: - // http://netlib.org/linalg/html_templates/node92.html - // or - // http://en.wikipedia.org/wiki/Sparse_matrix - // #Compressed_sparse_column_.28CSC_or_CCS.29 - // - // - // Templates: - // Index type for index variables - // Scalar type for floating point variables (gets cast to double?) - // Input: - // n number of variables, i.e. size of x - // Qi vector of qnnz row indices of non-zeros in LOWER TRIANGLE ONLY of - // Q⁰ - // Qj vector of qnnz column indices of non-zeros in LOWER TRIANGLE ONLY - // of Q⁰ - // Qv vector of qnnz values of non-zeros in LOWER TRIANGLE ONLY of Q⁰, - // such that: - // - // Q⁰(Qi[k],Qj[k]) = Qv[k] for k ∈ [0,Qnnz-1], where Qnnz is the - // - // number of non-zeros in Q⁰ - // c (optional) vector of n values of c, transpose of coefficient row - // vector of linear terms, EMPTY means c == 0 - // cf (ignored) value of constant term in objective, 0 means cf == 0, so - // optional only in the sense that it is mandatory - // m number of constraints, therefore also number of rows in linear - // constraint coefficient matrix A, and in linear constraint bound - // vectors lc and uc - // Av vector of non-zero values of A, in column compressed order - // Ari vector of row indices corresponding to non-zero values of A, - // Acp vector of indices into Ari and Av of the first entry for each - // column of A, size(Acp) = (# columns of A) + 1 = n + 1 - // lc vector of m linear constraint lower bounds - // uc vector of m linear constraint upper bounds - // lx vector of n constant lower bounds - // ux vector of n constant upper bounds - // Output: - // x vector of size n to hold output of optimization - // Return: - // true only if optimization was successful with no errors - // - // Note: All indices are 0-based - // - template <typename Index, typename Scalar> - IGL_INLINE bool mosek_quadprog( - const Index n, - /* mosek won't allow this to be const*/ std::vector<Index> & Qi, - /* mosek won't allow this to be const*/ std::vector<Index> & Qj, - /* mosek won't allow this to be const*/ std::vector<Scalar> & Qv, - const std::vector<Scalar> & c, - const Scalar cf, - const Index m, - /* mosek won't allow this to be const*/ std::vector<Scalar> & Av, - /* mosek won't allow this to be const*/ std::vector<Index> & Ari, - const std::vector<Index> & Acp, - const std::vector<Scalar> & lc, - const std::vector<Scalar> & uc, - const std::vector<Scalar> & lx, - const std::vector<Scalar> & ux, - MosekData & mosek_data, - std::vector<Scalar> & x); - // Wrapper with Eigen elements - // - // Inputs: - // Q n by n square quadratic coefficients matrix **only lower triangle - // is used**. - // c n-long vector of linear coefficients - // cf constant coefficient - // A m by n square linear coefficienst matrix of inequality constraints - // lc m-long vector of lower bounds for linear inequality constraints - // uc m-long vector of upper bounds for linear inequality constraints - // lx n-long vector of lower bounds - // ux n-long vector of upper bounds - // mosek_data parameters struct - // Outputs: - // x n-long solution vector - // Returns true only if optimization finishes without error - // - IGL_INLINE bool mosek_quadprog( - const Eigen::SparseMatrix<double> & Q, - const Eigen::VectorXd & c, - const double cf, - const Eigen::SparseMatrix<double> & A, - const Eigen::VectorXd & lc, - const Eigen::VectorXd & uc, - const Eigen::VectorXd & lx, - const Eigen::VectorXd & ux, - MosekData & mosek_data, - Eigen::VectorXd & x); - } -} - -#ifndef IGL_STATIC_LIBRARY -# include "mosek_quadprog.cpp" -#endif - -#endif |