diff --git a/.gitignore b/.gitignore
index dd20ea0d..77d6357a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,7 @@
+cmake-build-debug/
+cmake-build-release/
+build/
+.idea/
 examples/sampling/build/CMakeCache.txt
 examples/sampling/build/CMakeFiles/3.5.2/CMakeCCompiler.cmake
 examples/sampling/build/CMakeFiles/3.5.2/CMakeCXXCompiler.cmake
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 00000000..328c7e85
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,6 @@
+cmake_minimum_required(VERSION 3.20)
+project(Directional)
+set(CMAKE_CXX_STANDARD 17)
+
+add_library(${PROJECT_NAME} INTERFACE include)
+target_include_directories(${PROJECT_NAME} INTERFACE include)
diff --git a/include/directional/integrate.h b/include/directional/integrate.h
index 6017d63d..11f50aa3 100644
--- a/include/directional/integrate.h
+++ b/include/directional/integrate.h
@@ -127,153 +127,8 @@ namespace directional
         //reducing constraintMat
         SparseQR<SparseMatrix<double>, COLAMDOrdering<int> > qrsolver;
         SparseMatrix<double> Cfull = intData.constraintMat * intData.linRedMat * intData.singIntSpanMat * intData.intSpanMat;
-        if (Cfull.rows()!=0){
-            qrsolver.compute(Cfull.transpose());
-            int CRank = qrsolver.rank();
-
-            //creating sliced permutation matrix
-            VectorXi PIndices = qrsolver.colsPermutation().indices();
-
-            vector<Triplet<double> > CTriplets;
-            for(int k = 0; k < Cfull.outerSize(); ++k)
-            {
-                for(SparseMatrix<double>::InnerIterator it(Cfull, k); it; ++it)
-                {
-                    for(int j = 0; j < CRank; j++)
-                        if(it.row() == PIndices(j))
-                            CTriplets.emplace_back(j, it.col(), it.value());
-                }
-            }
-
-            Cfull.resize(CRank, Cfull.cols());
-            Cfull.setFromTriplets(CTriplets.begin(), CTriplets.end());
-        }
-        SparseMatrix<double> var2AllMat;
-        VectorXd fullx(numVars); fullx.setZero();
-        for(int intIter = 0; intIter < fixedMask.sum(); intIter++)
-        {
-            //the non-fixed variables to all variables
-            var2AllMat.resize(numVars, numVars - alreadyFixed.sum());
-            int varCounter = 0;
-            vector<Triplet<double> > var2AllTriplets;
-            for(int i = 0; i < numVars; i++)
-            {
-                if (!alreadyFixed(i)){
-                    //for (int j=0;j<intData.d;j++)
-                    var2AllTriplets.emplace_back(i, varCounter++, 1.0);
-                }
-
-            }
-            var2AllMat.setFromTriplets(var2AllTriplets.begin(), var2AllTriplets.end());
-
-            SparseMatrix<double> Epart = Efull * var2AllMat;
-            VectorXd torhs = -Efull * fixedValues;
-            SparseMatrix<double> EtE = Epart.transpose() * M1 * Epart;
-            SparseMatrix<double> Cpart = Cfull * var2AllMat;
-
-            //reducing rank on Cpart
-            int CpartRank=0;
-            VectorXi PIndices(0);
-            if (Cpart.rows()!=0){
-                qrsolver.compute(Cpart.transpose());
-                CpartRank = qrsolver.rank();
-
-                //creating sliced permutation matrix
-                PIndices = qrsolver.colsPermutation().indices();
-
-                vector<Triplet<double> > CPartTriplets;
-
-                for(int k = 0; k < Cpart.outerSize(); ++k)
-                {
-                    for (SparseMatrix<double>::InnerIterator it(Cpart, k); it; ++it)
-                    {
-                        for (int j = 0; j < CpartRank; j++)
-                            if (it.row() == PIndices(j))
-                                CPartTriplets.emplace_back(j, it.col(), it.value());
-                    }
-                }
-
-                Cpart.resize(CpartRank, Cpart.cols());
-                Cpart.setFromTriplets(CPartTriplets.begin(), CPartTriplets.end());
-            }
-            SparseMatrix<double> A(EtE.rows()+ Cpart.rows(), EtE.rows() + Cpart.rows());
-
-            vector<Triplet<double>> ATriplets;
-            for(int k = 0; k < EtE.outerSize(); ++k)
-            {
-                for (SparseMatrix<double>::InnerIterator it(EtE, k); it; ++it)
-                    ATriplets.push_back(Triplet<double>(it.row(), it.col(), it.value()));
-            }
-
-            for(int k = 0; k < Cpart.outerSize(); ++k)
-            {
-                for(SparseMatrix<double>::InnerIterator it(Cpart, k); it; ++it)
-                {
-                    ATriplets.emplace_back(it.row() + EtE.rows(), it.col(), it.value());
-                    ATriplets.emplace_back(it.col(), it.row() + EtE.rows(), it.value());
-                }
-            }
-
-            A.setFromTriplets(ATriplets.begin(), ATriplets.end());
-
-            //Right-hand side with fixed values
-            VectorXd b = VectorXd::Zero(EtE.rows() + Cpart.rows());
-            b.segment(0, EtE.rows())= Epart.transpose() * M1 * (gamma + torhs);
-            VectorXd bfull = -Cfull * fixedValues;
-            VectorXd bpart(CpartRank);
-            for(int k = 0; k < CpartRank; k++)
-                bpart(k)=bfull(PIndices(k));
-            b.segment(EtE.rows(), Cpart.rows()) = bpart;
-
-            SparseLU<SparseMatrix<double> > lusolver;
-            lusolver.compute(A);
-            if(lusolver.info() != Success){
-                if (intData.verbose)
-                    cout<<"LU decomposition failed!"<<endl;
-                return false;
-            }
-            x = lusolver.solve(b);
-
-            fullx = var2AllMat * x.head(numVars - alreadyFixed.sum()) + fixedValues;
-
-
-            if((alreadyFixed - fixedMask).sum() == 0)
-                break;
-
-            double minIntDiff = std::numeric_limits<double>::max();
-            int minIntDiffIndex = -1;
-            for (int i = 0; i < numVars; i++)
-            {
-                if ((fixedMask(i)) && (!alreadyFixed(i)))
-                {
-                    double currIntDiff =0;
-                    double func = fullx(i); //fullx.segment(intData.d*i,intData.d);
-                    //for (int j=0;j<intData.d;j++)
-                    currIntDiff += std::fabs(func - std::round(func));
-                    if (currIntDiff < minIntDiff)
-                    {
-                        minIntDiff = currIntDiff;
-                        minIntDiffIndex = i;
-                    }
-                }
-            }
-
-            if (minIntDiffIndex != -1)
-            {
-                alreadyFixed(minIntDiffIndex) = 1;
-                double func = fullx(minIntDiffIndex) ;
-                double funcInteger=std::round(func);
-                fixedValues(minIntDiffIndex) = /*pinvSymm*projMat**/funcInteger;
-            }
-
-            xprev.resize(x.rows() - 1);
-            varCounter = 0;
-            for(int i = 0; i < numVars; i++)
-                if (!alreadyFixed(i))
-                    xprev(varCounter++) = fullx(i);
-
-            xprev.tail(Cpart.rows()) = x.tail(Cpart.rows());
-        }
+        VectorXd fullx(numVars);
+        fullx.setZero();
 
         //the results are packets of N functions for each vertex, and need to be allocated for corners
         VectorXd NFunctionVec = intData.vertexTrans2CutMat * intData.linRedMat * intData.singIntSpanMat * intData.intSpanMat * fullx;