Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

removed previous iterative rounding codes #75

Closed
wants to merge 1 commit into from
Closed

removed previous iterative rounding codes #75

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Expand Down
6 changes: 6 additions & 0 deletions CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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)
149 changes: 2 additions & 147 deletions include/directional/integrate.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down