Some redundant stuff.

Former-commit-id: ee50eeb9262a4876fb07b04719cb5cd7e4929355

include/igl/DiscreteShellsTraits.h

@@ -1,315 +0,0 @@
-// This file is part of libhedra, a library for polyhedral mesh processing
-//
-// Copyright (C) 2016 Amir Vaxman <avaxman@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 HEDRA_DISCRETE_SHELLS_TRAITS_H
-#define HEDRA_DISCRETE_SHELLS_TRAITS_H
-#include <igl/igl_inline.h>
-#include <igl/harmonic.h>
-#include <Eigen/Core>
-#include <string>
-#include <vector>
-#include <cstdio>
-#include <set>
-
-
-namespace hedra { namespace optimization {
-        
-        //this class is a traits class for optimization of discrete shells deformation by given positional constraints. It is an implementation on [Froehlich and Botsch 2012] for general polyhedral meshes, using a triangulation of them
-    
-        //the solution vector is assumed to be arranged as xyzxyzxyz... where each triplet is a coordinate of the free vertices.
-    
-        class DiscreteShellsTraits{
-        public:
-            
-            //Requisites of the traits class
-            Eigen::VectorXi JRows, JCols;  //rows and column indices for the jacobian matrix
-            Eigen::VectorXd JVals;         //values for the jacobian matrix.
-            int xSize;                  //size of the solution
-            Eigen::VectorXd EVec;          //energy vector
-            
-            //These are for the the full Jacobian matrix without removing the handles
-            Eigen::VectorXi fullJRows, fullJCols;
-            Eigen::VectorXd fullJVals;
-            Eigen::MatrixXi flapVertexIndices;  //vertices (i,j,k,l) of a flap on edge e=(k,i) where the triangles are f=(i,j,k) and g=(i,k,l)
-            Eigen::MatrixXi EV;
-            Eigen::VectorXi h;              //list of handles
-            Eigen::MatrixXd qh;             //#h by 3 positions
-            Eigen::VectorXi a2x;            //from the entire set of variables "a" to the free variables in the optimization "x".
-            Eigen::VectorXi colMap;         //raw map of a2x into the columns from fullJCols into JCols
-            Eigen::VectorXd origLengths;    //the original edge lengths.
-            Eigen::VectorXd origDihedrals;  //Original dihedral angles
-            Eigen::MatrixXd VOrig;          //original positions
-            Eigen::MatrixXi T;              //triangles
-            Eigen::VectorXd Wd, Wl;         //geometric weights for the lengths and the dihedral angles.
-            double bendCoeff, lengthCoeff;  //coefficients for the different terms
-            
-            //Eigen::VectorXd x0;                 //the initial solution to the optimization
-            Eigen::MatrixXd fullSolution;       //The final solution of the last optimization
-            
-            void init(const Eigen::MatrixXd& _VOrig,
-                      const Eigen::MatrixXi& _T,
-                      const Eigen::VectorXi& _h,
-                      const Eigen::MatrixXi& _EV,
-                      const Eigen::MatrixXi& ET,
-                      const Eigen::MatrixXi& ETi,
-                      const Eigen::VectorXi& innerEdges){
-                
-                using namespace std;
-                using namespace Eigen;
-                
-                std::set<pair<int, int> > edgeIndicesList;
-                
-                VOrig=_VOrig;
-                T=_T;
-                h=_h;
-                EV=_EV;
-                lengthCoeff=10.0;
-                bendCoeff=1.0;
-                
-                //lengths of edges and diagonals
-                flapVertexIndices.resize(innerEdges.size(),4);
-                
-                //dihedral angles
-                for (int i=0;i<innerEdges.size();i++){
-                    int f=ET(innerEdges(i),0);
-                    int g=ET(innerEdges(i),1);
-                 
-                    int vi=EV(innerEdges(i), 1);
-                    int vj=T(f,(ETi(innerEdges(i),0)+2)%3);
-                    int vk=EV(innerEdges(i),0);
-                    int vl=T(g,(ETi(innerEdges(i),1)+2)%3);
-                    flapVertexIndices.row(i)<<vi,vj,vk,vl;
-                 }
-                
-                //across edges
-                EVec.resize(EV.rows()+flapVertexIndices.rows());
-                origLengths.resize(EV.rows());
-                origDihedrals.resize(flapVertexIndices.rows());
-                Wl.resize(EV.rows());
-                Wd.resize(flapVertexIndices.rows());
-                
-                
-                for (int i=0;i<EV.rows();i++){
-                    origLengths(i)=(VOrig.row(EV(i,1))-VOrig.row(EV(i,0))).norm();
-                    Wl(i)=1.0/origLengths(i);
-                }
-                
-                for (int i=0;i<flapVertexIndices.rows();i++){
-                    RowVector3d eji=VOrig.row(flapVertexIndices(i,0))-VOrig.row(flapVertexIndices(i,1));
-                    RowVector3d ejk=VOrig.row(flapVertexIndices(i,2))-VOrig.row(flapVertexIndices(i,1));
-                    RowVector3d eli=VOrig.row(flapVertexIndices(i,0))-VOrig.row(flapVertexIndices(i,3));
-                    RowVector3d elk=VOrig.row(flapVertexIndices(i,2))-VOrig.row(flapVertexIndices(i,3));
-                    RowVector3d eki=VOrig.row(flapVertexIndices(i,0))-VOrig.row(flapVertexIndices(i,2));
-                    
-                    RowVector3d n1 = (ejk.cross(eji));
-                    RowVector3d n2 = (eli.cross(elk));
-                    double sign=((n1.cross(n2)).dot(eki) >= 0 ? 1.0 : -1.0);
-                    double sinHalf=sign*sqrt((1.0-n1.normalized().dot(n2.normalized()))/2.0);
-                    origDihedrals(i)=2.0*asin(sinHalf);
-                    double areaSum=(n1.norm()+n2.norm())/2.0;
-                    Wd(i)=eki.norm()/sqrt(areaSum);
-                }
-                
-                
-                //creating the Jacobian pattern
-                xSize=3*(VOrig.rows()-h.size());
-                fullJRows.resize(6*EV.rows()+12*flapVertexIndices.rows());
-                fullJCols.resize(6*EV.rows()+12*flapVertexIndices.rows());
-                fullJVals.resize(6*EV.rows()+12*flapVertexIndices.rows());
-                
-                
-                //Jacobian indices for edge lengths
-                for (int i=0;i<EV.rows();i++){
-                    fullJRows.segment(6*i,6).setConstant(i);
-                    fullJCols.segment(6*i,3)<<3*EV(i,0),3*EV(i,0)+1,3*EV(i,0)+2;
-                    fullJCols.segment(6*i+3,3)<<3*EV(i,1),3*EV(i,1)+1,3*EV(i,1)+2;
-                }
-                
-                //Jacobian indices for dihedral angles
-                for (int i=0;i<flapVertexIndices.rows();i++){
-                    fullJRows.segment(6*EV.rows()+12*i,12).setConstant(EV.rows()+i);
-                    for (int k=0;k<4;k++)
-                        fullJCols.segment(6*EV.rows()+12*i+3*k,3)<<3*flapVertexIndices(i,k),3*flapVertexIndices(i,k)+1,3*flapVertexIndices(i,k)+2;
-                }
-                
-                
-                a2x=Eigen::VectorXi::Zero(VOrig.rows());
-                int CurrIndex=0;
-                for (int i=0;i<h.size();i++)
-                    a2x(h(i))=-1;
-                
-                for (int i=0;i<VOrig.rows();i++)
-                    if (a2x(i)!=-1)
-                        a2x(i)=CurrIndex++;
-                
-                colMap.resize(3*VOrig.rows());
-                for (int i=0;i<VOrig.rows();i++)
-                    if (a2x(i)!=-1)
-                        colMap.segment(3*i,3)<<3*a2x(i),3*a2x(i)+1,3*a2x(i)+2;
-                    else
-                        colMap.segment(3*i,3)<<-1,-1,-1;
-                
-                //setting up the Jacobian rows and columns
-                int actualGradCounter=0;
-                for (int i=0;i<fullJCols.size();i++){
-                    if (colMap(fullJCols(i))!=-1)  //not a removed variable
-                        actualGradCounter++;
-                }
-                
-                JRows.resize(actualGradCounter);
-                JCols.resize(actualGradCounter);
-                JVals.resize(actualGradCounter);
-                
-                
-                actualGradCounter=0;
-                for (int i=0;i<fullJCols.size();i++){
-                    if (colMap(fullJCols(i))!=-1){  //not a removed variable
-                        JRows(actualGradCounter)=fullJRows(i);
-                        JCols(actualGradCounter++)=colMap(fullJCols(i));
-                    }
-                }
-            }
-            
-            //provide the initial solution to the solver
-            void initial_solution(Eigen::VectorXd& x0){
-                //using biharmonic deformation fields
-                Eigen::MatrixXd origHandlePoses(qh.rows(),3);
-                for (int i=0;i<qh.rows();i++)
-                    origHandlePoses.row(i)=VOrig.row(h(i));
-                Eigen::MatrixXd D;
-                Eigen::MatrixXd D_bc = qh - origHandlePoses;
-                igl::harmonic(VOrig,T,h,D_bc,2,D);
-                Eigen::MatrixXd V0 = VOrig+D;
-                for (int i=0;i<VOrig.rows();i++)
-                    if (a2x(i)!=-1)
-                        x0.segment(3*a2x(i),3)<<V0.row(i).transpose();
-                
-            }
-            
-            void pre_iteration(const Eigen::VectorXd& prevx){}
-            bool post_iteration(const Eigen::VectorXd& x){return false;  /*never stop after an iteration*/}
-            
-            
-            //updating the energy vector for a given current solution
-            void update_energy(const Eigen::VectorXd& x){
-                
-                using namespace std;
-                using namespace Eigen;
-                
-                MatrixXd fullx(xSize+h.size(),3);
-                for (int i=0;i<a2x.size();i++)
-                    if (a2x(i)!=-1)
-                        fullx.row(i)<<x.segment(3*a2x(i),3).transpose();
-                
-                for (int i=0;i<h.size();i++)
-                    fullx.row(h(i))=qh.row(i);
-                
-                fullJVals.setZero();
-                for (int i=0;i<EV.rows();i++)
-                    EVec(i)=((fullx.row(EV(i,1))-fullx.row(EV(i,0))).norm()-origLengths(i))*Wl(i)*lengthCoeff;
-                
-                
-                for (int i=0;i<flapVertexIndices.rows();i++){
-                    RowVector3d eji=fullx.row(flapVertexIndices(i,0))-fullx.row(flapVertexIndices(i,1));
-                    RowVector3d ejk=fullx.row(flapVertexIndices(i,2))-fullx.row(flapVertexIndices(i,1));
-                    RowVector3d eli=fullx.row(flapVertexIndices(i,0))-fullx.row(flapVertexIndices(i,3));
-                    RowVector3d elk=fullx.row(flapVertexIndices(i,2))-fullx.row(flapVertexIndices(i,3));
-                    RowVector3d eki=fullx.row(flapVertexIndices(i,0))-fullx.row(flapVertexIndices(i,2));
-                    
-                    RowVector3d n1 = (ejk.cross(eji));
-                    RowVector3d n2 = (eli.cross(elk));
-                    double sign=((n1.cross(n2)).dot(eki) >= 0 ? 1.0 : -1.0);
-                    double dotn1n2=1.0-n1.normalized().dot(n2.normalized());
-                    if (dotn1n2<0.0) dotn1n2=0.0;  //sanitizing
-                    double sinHalf=sign*sqrt(dotn1n2/2.0);
-                    if (sinHalf>1.0) sinHalf=1.0; if (sinHalf<-1.0) sinHalf=-1.0;
-                    EVec(EV.rows()+i)=(2.0*asin(sinHalf)-origDihedrals(i))*Wd(i)*bendCoeff;
-                }
-                
-                for (int i=0;i<EVec.size();i++)
-                    if (isnan(EVec(i)))
-                        cout<<"nan in EVec("<<i<<")"<<endl;
-            }
-            
-            
-            //update the jacobian values for a given current solution
-            void update_jacobian(const Eigen::VectorXd& x){
-                using namespace std;
-                using namespace Eigen;
-                
-                MatrixXd fullx(xSize+h.size(),3);
-                for (int i=0;i<a2x.size();i++)
-                    if (a2x(i)!=-1)
-                        fullx.row(i)<<x.segment(3*a2x(i),3).transpose();
-                
-                for (int i=0;i<h.size();i++)
-                    fullx.row(h(i))=qh.row(i);
-                
-                fullJVals.setZero();
-                for (int i=0;i<EV.rows();i++){
-                   
-                    RowVector3d normedEdgeVector=(fullx.row(EV(i,1))-fullx.row(EV(i,0))).normalized();
-                    fullJVals.segment(6*i,3)<<-normedEdgeVector.transpose()*Wl(i)*lengthCoeff;
-                    fullJVals.segment(6*i+3,3)<<normedEdgeVector.transpose()*Wl(i)*lengthCoeff;
-                }
-                
-                for (int i=0;i<flapVertexIndices.rows();i++){
-                    RowVector3d eji=fullx.row(flapVertexIndices(i,0))-fullx.row(flapVertexIndices(i,1));
-                    RowVector3d ejk=fullx.row(flapVertexIndices(i,2))-fullx.row(flapVertexIndices(i,1));
-                    RowVector3d eli=fullx.row(flapVertexIndices(i,0))-fullx.row(flapVertexIndices(i,3));
-                    RowVector3d elk=fullx.row(flapVertexIndices(i,2))-fullx.row(flapVertexIndices(i,3));
-                    RowVector3d eki=fullx.row(flapVertexIndices(i,0))-fullx.row(flapVertexIndices(i,2));
-                    
-                    RowVector3d n1 = (ejk.cross(eji));
-                    RowVector3d n2 = (eli.cross(elk));
-                    double sign=((n1.cross(n2)).dot(eki) >= 0 ? 1.0 : -1.0);
-                    double dotn1n2=1.0-n1.normalized().dot(n2.normalized());
-                    if (dotn1n2<0.0) dotn1n2=0.0;  //sanitizing
-                    double sinHalf=sign*sqrt(dotn1n2/2.0);
-                    if (sinHalf>1.0) sinHalf=1.0; if (sinHalf<-1.0) sinHalf=-1.0;
-                    
-                    fullJVals.segment(6*EV.rows()+12*i,3)<<(Wd(i)*((ejk.dot(-eki)/(n1.squaredNorm()*eki.norm()))*n1+(elk.dot(-eki)/(n2.squaredNorm()*eki.norm()))*n2)).transpose()*bendCoeff;
-                    fullJVals.segment(6*EV.rows()+12*i+3,3)<<(Wd(i)*(-eki.norm()/n1.squaredNorm())*n1).transpose()*bendCoeff;
-                    fullJVals.segment(6*EV.rows()+12*i+6,3)<<(Wd(i)*((eji.dot(eki)/(n1.squaredNorm()*eki.norm()))*n1+(eli.dot(eki)/(n2.squaredNorm()*eki.norm()))*n2)).transpose()*bendCoeff;
-                    fullJVals.segment(6*EV.rows()+12*i+9,3)<<(Wd(i)*(-eki.norm()/n2.squaredNorm())*n2).transpose()*bendCoeff;
-                    
-                }
-                
-
-                int actualGradCounter=0;
-                for (int i=0;i<fullJCols.size();i++)
-                    if (colMap(fullJCols(i))!=-1)  //not a removed variable
-                        JVals(actualGradCounter++)=fullJVals(i);
-                
-                for (int i=0;i<JVals.size();i++)
-                    if (isnan(JVals(i)))
-                        cout<<"nan in JVals("<<i<<")"<<endl;
-            }
-
-            
-            
-            bool post_optimization(const Eigen::VectorXd& x){
-                fullSolution.conservativeResize(a2x.size(),3);
-                for (int i=0;i<a2x.size();i++)
-                    if (a2x(i)!=-1)
-                        fullSolution.row(i)<<x.segment(3*a2x(i),3).transpose();
-                
-                for (int i=0;i<h.size();i++)
-                    fullSolution.row(h(i))=qh.row(i);
-                
-                return true;  //stop optimization after this
-            }
-            
-            DiscreteShellsTraits(){}
-            ~DiscreteShellsTraits(){}
-        };
-        
-        
-    } }
-
-
-#endif

include/igl/LMSolver.h

@@ -1,350 +0,0 @@
-// This file is part of libhedra, a library for polyhedral mesh processing
-//
-// Copyright (C) 2016 Amir Vaxman <avaxman@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 HEDRA_LEVENBERG_MARQUADT_SOLVER_H
-#define HEDRA_LEVENBERG_MARQUADT_SOLVER_H
-#include <igl/igl_inline.h>
-#include <igl/sortrows.h>
-#include <igl/speye.h>
-#include <Eigen/Core>
-#include <string>
-#include <vector>
-#include <cstdio>
-#include <iostream>
-
-namespace hedra {
-    namespace optimization
-    {
-        
-        template<class LinearSolver, class SolverTraits>
-        class LMSolver{
-        public:
-            Eigen::VectorXd x;      //current solution; always updated
-            Eigen::VectorXd prevx;  //the solution of the previous iteration
-            Eigen::VectorXd x0;     //the initial solution to the system
-            Eigen::VectorXd d;             //the direction taken.
-            Eigen::VectorXd currEnergy;    //the current value of the energy
-            Eigen::VectorXd prevEnergy;    //the previous value of the energy
-            
-            Eigen::VectorXi HRows, HCols;  //(row,col) pairs for H=J^T*J matrix
-            Eigen::VectorXd HVals;      //values for H matrix
-            Eigen::MatrixXi S2D;        //single J to J^J indices
-
-            LinearSolver* LS;
-            SolverTraits* ST;
-            int maxIterations;
-            double xTolerance;
-            double fooTolerance;
-            
-            /*void TestMatrixOperations(){
-             
-                using namespace Eigen;
-                using namespace std;
-                int RowSize=1000;
-                int ColSize=1000;
-                int TriSize=4000;
-             
-                VectorXi Rows;
-                VectorXi Cols;
-                VectorXd Vals=VectorXd::Random(TriSize);
-             
-                VectorXd dRows=VectorXd::Random(TriSize);
-                VectorXd dCols=VectorXd::Random(TriSize);
-             
-                cout<<"dRows Range: "<<dRows.minCoeff()<<","<<dRows.maxCoeff()<<endl;
-             
-                Rows=((dRows+MatrixXd::Constant(dRows.size(),1,1.0))*500.0).cast<int>();
-                Cols=((dCols+MatrixXd::Constant(dRows.size(),1,1.0))*500.0).cast<int>();
-                VectorXi SortRows, stub;
-                VectorXi MRows, MCols;
-                VectorXd MVals;
-             
-             igl::sortrows(Rows, true, SortRows, stub);
-             
-             MatrixXi S2D;
-                
-            double miu=15.0;
-             
-             MatrixPattern(SortRows, Cols,  MRows, MCols, S2D);
-             MVals.resize(MRows.size());
-             MatrixValues(MRows, MCols, Vals, S2D, miu, MVals);
- 
-             //testing multiplyadjoint
-             SparseMatrix<double> M(SortRows.maxCoeff()+1,Cols.maxCoeff()+1);
-             
-             //cout<<"Our I"<<I<<endl;
-             //cout<<"Our J"<<J<<endl;
-             //cout<<"Our S"<<S<<endl;
-             
-             
-             vector<Triplet<double> > Tris;
-             
-             for (int i=0;i<SortRows.size();i++)
-             Tris.push_back(Triplet<double>(SortRows(i), Cols(i), Vals(i)));
-             
-             M.setFromTriplets(Tris.begin(), Tris.end());
-             Tris.clear();
-                
-                
-                SparseMatrix<double> I;
-                igl::speye(M.cols(),M.cols(),I);
-             
-             SparseMatrix<double> MtM1=M.transpose()*M+miu*I;
-             SparseMatrix<double> MtM2(MtM1.rows(), MtM1.cols());
-             for (int i=0;i<MRows.size();i++)
-                 Tris.push_back(Triplet<double>(MRows(i), MCols(i), MVals(i)));
-             
-             MtM2.setFromTriplets(Tris.begin(), Tris.end());
-             
-             bool Discrepancy=false;
-             SparseMatrix<double> DiffMat=MtM1-MtM2;
-             for (int k=0; k<DiffMat.outerSize(); ++k){
-             for (SparseMatrix<double>::InnerIterator it(DiffMat,k); it; ++it){
-             if ((abs(it.value())>10e-6)&&(it.row()<=it.col())){
-             cout<<"Discrepancy at ("<<it.row()<<","<<it.col()<<"): "<<it.value()<<endl;
-             cout<<"MtM Our Values:"<<MtM2.coeffRef(it.row(),it.col())<<endl;
-             cout<<"MtM Evaluated:"<<MtM1.coeffRef(it.row(),it.col())<<endl;
-             Discrepancy=true;
-             }
-             }
-             }
-             if (!Discrepancy) cout<<"Matrices are equal!"<<endl;
-             }*/
-
-            
-            //Input: pattern of matrix M by (iI,iJ) representation
-            //Output: pattern of matrix M^T*M by (oI, oJ) representation
-            //        map between values in the input to values in the output (Single2Double). The map is aggregating values from future iS to oS
-            //prerequisite: iI are sorted by rows (not necessary columns)
-            void MatrixPattern(const Eigen::VectorXi& iI,
-                               const Eigen::VectorXi& iJ,
-                               Eigen::VectorXi& oI,
-                               Eigen::VectorXi& oJ,
-                               Eigen::MatrixXi& S2D)
-            {
-                int CurrTri=0;
-                using namespace Eigen;
-                std::vector<int> oIlist;
-                std::vector<int> oJlist;
-                std::vector<std::pair<int, int> > S2Dlist;
-                do{
-                    int CurrRow=iI(CurrTri);
-                    int NumCurrTris=0;
-                    while ((CurrTri+NumCurrTris<iI.size())&&(iI(CurrTri+NumCurrTris)==CurrRow))
-                        NumCurrTris++;
-                    
-                    for (int i=CurrTri;i<CurrTri+NumCurrTris;i++){
-                        for (int j=CurrTri;j<CurrTri+NumCurrTris;j++){
-                            if (iJ(j)>=iJ(i)){
-                                oIlist.push_back(iJ(i));
-                                oJlist.push_back(iJ(j));
-                                S2Dlist.push_back(std::pair<int,int>(i,j));
-                            }
-                        }
-                    }
-                    CurrTri+=NumCurrTris;
-                }while (CurrTri!=iI.size());
-                
-                
-                //triplets for miu
-                for (int i=0;i<iJ.maxCoeff()+1;i++){
-                    oIlist.push_back(i);
-                    oJlist.push_back(i);
-                }
-                
-                oI.resize(oIlist.size());
-                oJ.resize(oJlist.size());
-                S2D.resize(S2Dlist.size(),2);
-                
-                for (int i=0;i<oIlist.size();i++){
-                    oI(i)=oIlist[i];
-                    oJ(i)=oJlist[i];
-                }
-                for (int i=0;i<S2Dlist.size();i++)
-                    S2D.row(i)<<S2Dlist[i].first, S2Dlist[i].second;
-                
-            }
-            
-            //returns the values of M^T*M+miu*I by multiplication and aggregating from Single2double list.
-            //prerequisite - oS is allocated
-            void MatrixValues(const Eigen::VectorXi& oI,
-                              const Eigen::VectorXi& oJ,
-                              const Eigen::VectorXd& iS,
-                              const Eigen::MatrixXi& S2D,
-                              const double miu,
-                              Eigen::VectorXd& oS)
-            {
-                for (int i=0;i<S2D.rows();i++)
-                    oS(i)=iS(S2D(i,0))*iS(S2D(i,1));
-                
-                //adding miu*I
-                for (int i=S2D.rows();i<oI.rows();i++)
-                    oS(i)=miu;
-                
-            }
-            
-            //returns M^t*ivec by (I,J,S) representation
-            void MultiplyAdjointVector(const Eigen::VectorXi& iI,
-                                       const Eigen::VectorXi& iJ,
-                                       const Eigen::VectorXd& iS,
-                                       const Eigen::VectorXd& iVec,
-                                       Eigen::VectorXd& oVec)
-            {
-                oVec.setZero();
-                for (int i=0;i<iI.size();i++)
-                    oVec(iJ(i))+=iS(i)*iVec(iI(i));
-            }
-            
-            
-        public:
-            
-            LMSolver(){};
-            
-            void init(LinearSolver* _LS,
-                      SolverTraits* _ST,
-                      int _maxIterations=100,
-                      double _xTolerance=10e-9,
-                      double _fooTolerance=10e-9){
-                
-                LS=_LS;
-                ST=_ST;
-                maxIterations=_maxIterations;
-                xTolerance=_xTolerance;
-                fooTolerance=_fooTolerance;
-                //analysing pattern
-                MatrixPattern(ST->JRows, ST->JCols,HRows,HCols,S2D);
-                HVals.resize(HRows.size());
-                
-                LS->analyze(HRows,HCols, true);
-                
-                d.resize(ST->xSize);
-                x.resize(ST->xSize);
-                x0.resize(ST->xSize);
-                prevx.resize(ST->xSize);
-                currEnergy.resize(ST->EVec.size());
-                prevEnergy.resize(ST->EVec.size());
-                
-                //TestMatrixOperations();
-            }
-            
-            
-            bool solve(const bool verbose) {
-                
-                using namespace Eigen;
-                using namespace std;
-                ST->initial_solution(x0);
-                prevx<<x0;
-                int currIter=0;
-                bool stop=false;
-                double currError, prevError;
-                VectorXd rhs(ST->xSize);
-                VectorXd direction;
-                if (verbose)
-                    cout<<"******Beginning Optimization******"<<endl;
-                
-                double tau=10e-3;
-                
-                //estimating initial miu
-                double miu=0.0;
-                ST->update_jacobian(prevx);
-                MatrixValues(HRows, HCols, ST->JVals, S2D, miu, HVals);
-                for (int i=0;i<HRows.size();i++)
-                    if (HRows(i)==HCols(i))  //on the diagonal
-                        miu=(miu < HVals(i) ? HVals(i) : miu);
-                miu*=tau;
-                double initmiu=miu;
-                cout<<"initial miu: "<<miu<<endl;
-                double beta=2.0;
-                double nu=beta;
-                double gamma=3.0;
-                do{
-                    currIter=0;
-                    stop=false;
-                    do{
-                        ST->pre_iteration(prevx);
-                        ST->update_energy(prevx);
-                        ST->update_jacobian(prevx);
-                        if (verbose)
-                            cout<<"Initial Energy for Iteration "<<currIter<<": "<<ST->EVec.template squaredNorm()<<endl;
-                        MatrixValues(HRows, HCols, ST->JVals, S2D,  miu, HVals);
-                        MultiplyAdjointVector(ST->JRows, ST->JCols, ST->JVals, -ST->EVec, rhs);
-                        
-                        double firstOrderOptimality=rhs.template lpNorm<Infinity>();
-                        if (verbose)
-                            cout<<"firstOrderOptimality: "<<firstOrderOptimality<<endl;
-                        
-                        if (firstOrderOptimality<fooTolerance){
-                            x=prevx;
-                            if (verbose){
-                                cout<<"First-order optimality has been reached"<<endl;
-                                break;
-                            }
-                        }
-                        
-                        //solving to get the GN direction
-                        if(!LS->factorize(HVals, true)) {
-                            // decomposition failed
-                            cout<<"Solver Failed to factorize! "<<endl;
-                            return false;
-                        }
-                        
-                        MatrixXd mRhs=rhs;
-                        MatrixXd mDirection;
-                        LS->solve(mRhs,mDirection);
-                        direction=mDirection.col(0);
-                        if (verbose)
-                            cout<<"direction magnitude: "<<direction.norm()<<endl;
-                        if (direction.norm() < xTolerance * prevx.norm()){
-                            x=prevx;
-                            if (verbose)
-                                cout<<"Stopping since direction magnitude small."<<endl;
-                            break;
-                        }
-                        VectorXd tryx=prevx+direction;
-                        ST->update_energy(prevx);
-                        double prevE=ST->EVec.squaredNorm();
-                        ST->update_energy(tryx);
-                        double currE=ST->EVec.squaredNorm();
-                        
-                        double rho=(prevE-currE)/(direction.dot(miu*direction+rhs));
-                        if (rho>0){
-                            x=tryx;
-                            //if (verbose){
-                                //cout<<"Energy: "<<currE<<endl;
-                            //    cout<<"1.0-(beta-1.0)*pow(2.0*rho-1.0,3): "<<1.0-(beta-1.0)*pow(2.0*rho-1.0,3)<<endl;
-                            //}
-                            miu*=(1.0/gamma > 1.0-(beta-1.0)*pow(2.0*rho-1.0,3) ? 1.0/gamma : 1.0-(beta-1.0)*pow(2.0*rho-1.0,3));
-                            nu=beta;
-                            //if (verbose)
-                            //    cout<<"rho, miu, nu: "<<rho<<","<<miu<<","<<nu<<endl;
-                        } else {
-                            x=prevx;
-                            miu = miu*nu;
-                            nu=2*nu;
-                            //if (verbose)
-                            //    cout<<"rho, miu, nu: "<<rho<<","<<miu<<","<<nu<<endl;
-                        }
-                                               
-                        //The SolverTraits can order the optimization to stop by giving "true" of to continue by giving "false"
-                        if (ST->post_iteration(x)){
-                            if (verbose)
-                                cout<<"ST->Post_iteration() gave a stop"<<endl;
-                            break;
-                        }
-                        currIter++;
-                        prevx=x;
-                    }while (currIter<=maxIterations);
-                }while (!ST->post_optimization(x));
-                return true;
-            }
-        };
-        
-    }
-}
-
-
-#endif

include/igl/SLSolver.h

@@ -1,184 +0,0 @@
-// This file is part of libhedra, a library for polyhedral mesh processing
-//
-// Copyright (C) 2016 Amir Vaxman <avaxman@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 HEDRA_SINGLE_LINEAR_SOLVER_H
-#define HEDRA_SINGLE_LINEAR_SOLVER_H
-#include <igl/igl_inline.h>
-#include <igl/sortrows.h>
-#include <igl/speye.h>
-#include <Eigen/Core>
-#include <string>
-#include <vector>
-#include <cstdio>
-#include <iostream>
-
-namespace hedra {
-    namespace optimization
-    {
-        
-        template<class LinearSolver, class SolverTraits>
-        class SLSolver{
-        public:
-            Eigen::VectorXd x;      //current solution; always updated
-            Eigen::VectorXd prevx;  //the solution of the previous iteration
-            Eigen::VectorXd x0;     //the initial solution to the system
-            Eigen::VectorXd d;             //the direction taken.
-            Eigen::VectorXd currEnergy;    //the current value of the energy
-            Eigen::VectorXd prevEnergy;    //the previous value of the energy
-            
-            Eigen::VectorXi HRows, HCols;  //(row,col) pairs for H=J^T*J matrix
-            Eigen::VectorXd HVals;      //values for H matrix
-            Eigen::MatrixXi S2D;        //single J to J^J indices
-
-            LinearSolver* LS;
-            SolverTraits* ST;
-            
-
-            //Input: pattern of matrix M by (iI,iJ) representation
-            //Output: pattern of matrix M^T*M by (oI, oJ) representation
-            //        map between values in the input to values in the output (Single2Double). The map is aggregating values from future iS to oS
-            //prerequisite: iI are sorted by rows (not necessary columns)
-            void MatrixPattern(const Eigen::VectorXi& iI,
-                               const Eigen::VectorXi& iJ,
-                               Eigen::VectorXi& oI,
-                               Eigen::VectorXi& oJ,
-                               Eigen::MatrixXi& S2D)
-            {
-                int CurrTri=0;
-                using namespace Eigen;
-                std::vector<int> oIlist;
-                std::vector<int> oJlist;
-                std::vector<std::pair<int, int> > S2Dlist;
-                do{
-                    int CurrRow=iI(CurrTri);
-                    int NumCurrTris=0;
-                    while ((CurrTri+NumCurrTris<iI.size())&&(iI(CurrTri+NumCurrTris)==CurrRow))
-                        NumCurrTris++;
-                    
-                    for (int i=CurrTri;i<CurrTri+NumCurrTris;i++){
-                        for (int j=CurrTri;j<CurrTri+NumCurrTris;j++){
-                            if (iJ(j)>=iJ(i)){
-                                oIlist.push_back(iJ(i));
-                                oJlist.push_back(iJ(j));
-                                S2Dlist.push_back(std::pair<int,int>(i,j));
-                            }
-                        }
-                    }
-                    CurrTri+=NumCurrTris;
-                }while (CurrTri!=iI.size());
-            
-                
-                oI.resize(oIlist.size());
-                oJ.resize(oJlist.size());
-                S2D.resize(S2Dlist.size(),2);
-                
-                for (int i=0;i<oIlist.size();i++){
-                    oI(i)=oIlist[i];
-                    oJ(i)=oJlist[i];
-                }
-                for (int i=0;i<S2Dlist.size();i++)
-                    S2D.row(i)<<S2Dlist[i].first, S2Dlist[i].second;
-                
-            }
-            
-            //returns the values of M^T*M+miu*I by multiplication and aggregating from Single2double list.
-            //prerequisite - oS is allocated
-            void MatrixValues(const Eigen::VectorXi& oI,
-                              const Eigen::VectorXi& oJ,
-                              const Eigen::VectorXd& iS,
-                              const Eigen::MatrixXi& S2D,
-                              Eigen::VectorXd& oS)
-            {
-                for (int i=0;i<S2D.rows();i++)
-                    oS(i)=iS(S2D(i,0))*iS(S2D(i,1));
-                
-            }
-            
-            //returns M^t*ivec by (I,J,S) representation
-            void MultiplyAdjointVector(const Eigen::VectorXi& iI,
-                                       const Eigen::VectorXi& iJ,
-                                       const Eigen::VectorXd& iS,
-                                       const Eigen::VectorXd& iVec,
-                                       Eigen::VectorXd& oVec)
-            {
-                oVec.setZero();
-                for (int i=0;i<iI.size();i++)
-                    oVec(iJ(i))+=iS(i)*iVec(iI(i));
-            }
-            
-            
-        public:
-            
-            SLSolver(){};
-            
-            void init(LinearSolver* _LS,
-                      SolverTraits* _ST){
-                
-                LS=_LS;
-                ST=_ST;
-                //analysing pattern
-                MatrixPattern(ST->JRows, ST->JCols,HRows,HCols,S2D);
-                HVals.resize(HRows.size());
-                
-                LS->analyze(HRows,HCols);
-                
-                d.resize(ST->xSize);
-                x.resize(ST->xSize);
-                x0.resize(ST->xSize);
-                prevx.resize(ST->xSize);
-                currEnergy.resize(ST->EVec.size());
-                prevEnergy.resize(ST->EVec.size());
-                
-                //TestMatrixOperations();
-            }
-            
-            
-            bool solve(const bool verbose) {
-                
-                using namespace Eigen;
-                using namespace std;
-                ST->initial_solution(x0);
-                prevx<<x0;
-                int currIter=0;
-                bool stop=false;
-                double currError, prevError;
-                VectorXd rhs(ST->xSize);
-                VectorXd direction;
-                if (verbose)
-                    cout<<"******Beginning Optimization******"<<endl;
-
-                ST->update_jacobian(prevx);
-                do{
-                    ST->update_energy(prevx);
-                    ST->update_jacobian(prevx);
-                    ST->pre_iteration(prevx);
-                    MatrixValues(HRows, HCols, ST->JVals, S2D, HVals);
-                    MultiplyAdjointVector(ST->JRows, ST->JCols, ST->JVals, -ST->EVec, rhs);
-                    
-                    //solving to get the GN direction
-                    if(!LS->factorize(HVals)) {
-                        // decomposition failed
-                        cout<<"Solver Failed to factorize! "<<endl;
-                        return false;
-                    }
-                        
-                    LS->solve(rhs,direction);
-                    x=prevx+direction;
-                    ST->update_energy(x);
-                    ST->update_jacobian(x);
-                    ST->post_iteration(x);
-                    prevx=x;
-                }while (!ST->post_optimization(x));
-                return true;
-            }
-        };
-        
-    }
-}
-
-
-#endif

include/igl/check_traits.h

@@ -1,88 +0,0 @@
-// This file is part of libhedra, a library for polyhedral mesh processing
-//
-// Copyright (C) 2016 Amir Vaxman <avaxman@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 HEDRA_CHECK_TRAITS_H
-#define HEDRA_CHECK_TRAITS_H
-#include <igl/igl_inline.h>
-#include <Eigen/Core>
-#include <string>
-#include <vector>
-#include <cstdio>
-
-namespace hedra {
-    namespace optimization
-    {
-        //This function checks the Jacobian of a traits class that is put for optimization by approximate finite differences, and reports the difference. It is important to use after coding, but it is not for the actual optimization process, since it is really brute-force and slow.
-        template<class SolverTraits>
-        void check_traits(SolverTraits& Traits){
-            using namespace Eigen;
-            using namespace std;
-            cout<<"WARNING: FE gradient checking, reducing performance!!!"<<endl;
-            cout<<"******************************************************"<<endl;
-            VectorXd CurrSolution(Traits.xSize);
-            CurrSolution<<VectorXd::Random(Traits.xSize, 1);
-            
-            cout<<"Solution Size: "<<Traits.xSize<<endl;
-            
-            Traits.update_energy(CurrSolution);
-            Traits.update_jacobian(CurrSolution);
-         
-            int MaxRow=Traits.JRows.maxCoeff()+1;
-            vector<Triplet<double> > GradTris;
-            
-            for (int i=0;i<Traits.JRows.size();i++)
-                GradTris.push_back(Triplet<double>(Traits.JRows(i), Traits.JCols(i), Traits.JVals(i)));
-            
-            
-            SparseMatrix<double> TraitGradient(MaxRow, CurrSolution.size());
-            TraitGradient.setFromTriplets(GradTris.begin(),GradTris.end());
-            
-            SparseMatrix<double> FEGradient(MaxRow, CurrSolution.size());
-            vector<Triplet<double> > FEGradientTris;
-            for (int i=0;i<CurrSolution.size();i++){
-                VectorXd vh(CurrSolution.size()); vh.setZero(); vh(i)=10e-5;
-                Traits.update_energy(CurrSolution+vh);
-                Traits.update_jacobian(CurrSolution+vh);
-                VectorXd EnergyPlus=Traits.EVec;
-                Traits.update_energy(CurrSolution-vh);
-                Traits.update_jacobian(CurrSolution-vh);
-                VectorXd EnergyMinus=Traits.EVec;
-                VectorXd CurrGradient=(EnergyPlus-EnergyMinus)/(2*10e-5);
-                //cout<<CurrGradient<<endl;
-                for (int j=0;j<CurrGradient.size();j++)
-                    if (abs(CurrGradient(j))>10e-7)
-                        FEGradientTris.push_back(Triplet<double>(j,i,CurrGradient(j)));
-            }
-            
-            FEGradient.setFromTriplets(FEGradientTris.begin(), FEGradientTris.end());
-            SparseMatrix<double> DiffMat=FEGradient-TraitGradient;
-            double maxcoeff=-32767.0;
-            int Maxi,Maxj;
-            for (int k=0; k<DiffMat.outerSize(); ++k)
-                for (SparseMatrix<double>::InnerIterator it(DiffMat,k); it; ++it){
-                    if (maxcoeff<abs(it.value())){
-                        maxcoeff=abs(it.value());
-                        Maxi=it.row();
-                        Maxj=it.col();
-                        
-                    }
-                    if (abs(it.value())>10e-7){
-                        cout<<"Gradient Discrepancy at: ("<<it.row()<<","<<it.col()<<") of "<<it.value()<<endl;
-                        cout<<"Our Value: "<<TraitGradient.coeffRef(it.row(), it.col())<<endl;
-                        cout<<"Computed Value: "<<FEGradient.coeffRef(it.row(), it.col())<<endl;
-                    }
-                }
-            cout<<"Maximum gradient difference: "<<maxcoeff<<endl;
-            cout<<"At Location: ("<<Maxi<<","<<Maxj<<")"<<endl;
-            
-            
-        }
-    }
-}
-
-
-#endif