disable printing
parent
f88c928ca0
commit
befe397f7a
|
|
@ -657,7 +657,6 @@ std::pair<boost::shared_ptr<GaussianConditional>,
|
||||||
// Combine and sort variable blocks in elimination order
|
// Combine and sort variable blocks in elimination order
|
||||||
JacobianFactor::shared_ptr jointFactor;
|
JacobianFactor::shared_ptr jointFactor;
|
||||||
try {
|
try {
|
||||||
cout << "JacobianFactor make_shared" << endl;
|
|
||||||
jointFactor = boost::make_shared<JacobianFactor>(factors, keys);
|
jointFactor = boost::make_shared<JacobianFactor>(factors, keys);
|
||||||
} catch (std::invalid_argument&) {
|
} catch (std::invalid_argument&) {
|
||||||
throw InvalidDenseElimination(
|
throw InvalidDenseElimination(
|
||||||
|
|
@ -665,8 +664,6 @@ std::pair<boost::shared_ptr<GaussianConditional>,
|
||||||
"involved in the provided factors.");
|
"involved in the provided factors.");
|
||||||
}
|
}
|
||||||
|
|
||||||
jointFactor->print("JointFactor0:");
|
|
||||||
|
|
||||||
// Do dense elimination
|
// Do dense elimination
|
||||||
if (jointFactor->model_)
|
if (jointFactor->model_)
|
||||||
jointFactor->model_ = jointFactor->model_->QR(jointFactor->Ab_.matrix());
|
jointFactor->model_ = jointFactor->model_->QR(jointFactor->Ab_.matrix());
|
||||||
|
|
@ -676,13 +673,9 @@ std::pair<boost::shared_ptr<GaussianConditional>,
|
||||||
// Zero below the diagonal
|
// Zero below the diagonal
|
||||||
jointFactor->Ab_.matrix().triangularView<Eigen::StrictlyLower>().setZero();
|
jointFactor->Ab_.matrix().triangularView<Eigen::StrictlyLower>().setZero();
|
||||||
|
|
||||||
factors.print("Factors to eliminate: ");
|
|
||||||
jointFactor->print("JointFactor1:");
|
|
||||||
|
|
||||||
// Split elimination result into conditional and remaining factor
|
// Split elimination result into conditional and remaining factor
|
||||||
GaussianConditional::shared_ptr conditional = jointFactor->splitConditional(
|
GaussianConditional::shared_ptr conditional = jointFactor->splitConditional(
|
||||||
keys.size());
|
keys.size());
|
||||||
cout << "JacobianFactor split conditoinal ok!" << endl;
|
|
||||||
|
|
||||||
return make_pair(conditional, jointFactor);
|
return make_pair(conditional, jointFactor);
|
||||||
}
|
}
|
||||||
|
|
|
||||||
|
|
@ -415,7 +415,7 @@ Constrained::shared_ptr Constrained::unit(size_t augmentedDim) const {
|
||||||
// It is Gram-Schmidt orthogonalization rather than Householder
|
// It is Gram-Schmidt orthogonalization rather than Householder
|
||||||
// Previously Diagonal::QR
|
// Previously Diagonal::QR
|
||||||
SharedDiagonal Constrained::QR(Matrix& Ab) const {
|
SharedDiagonal Constrained::QR(Matrix& Ab) const {
|
||||||
bool verbose = true;
|
bool verbose = false;
|
||||||
if (verbose) cout << "\nStarting Constrained::QR" << endl;
|
if (verbose) cout << "\nStarting Constrained::QR" << endl;
|
||||||
|
|
||||||
// get size(A) and maxRank
|
// get size(A) and maxRank
|
||||||
|
|
@ -428,34 +428,29 @@ SharedDiagonal Constrained::QR(Matrix& Ab) const {
|
||||||
|
|
||||||
Vector pseudo(m); // allocate storage for pseudo-inverse
|
Vector pseudo(m); // allocate storage for pseudo-inverse
|
||||||
Vector invsigmas = reciprocal(sigmas_);
|
Vector invsigmas = reciprocal(sigmas_);
|
||||||
gtsam::print(invsigmas, "invsigmas: ");
|
|
||||||
Vector weights = emul(invsigmas,invsigmas); // calculate weights once
|
Vector weights = emul(invsigmas,invsigmas); // calculate weights once
|
||||||
|
|
||||||
// Denote weights of inequality constraints with the negative sign
|
// Denote weights of inequality constraints with the negative sign
|
||||||
// TODO: might be slow!
|
// TODO: might be slow!
|
||||||
for (size_t s = 0; s<sigmas_.size(); ++s)
|
for (size_t s = 0; s<sigmas_.size(); ++s)
|
||||||
weights[s] = (sigmas_[s]<0)?-weights[s]:weights[s];
|
weights[s] = (sigmas_[s]<0)?-weights[s]:weights[s];
|
||||||
gtsam::print(weights, "weights: ");
|
if (verbose) gtsam::print(weights, "weights: ");
|
||||||
|
|
||||||
// We loop over all columns, because the columns that can be eliminated
|
// We loop over all columns, because the columns that can be eliminated
|
||||||
// are not necessarily contiguous. For each one, estimate the corresponding
|
// are not necessarily contiguous. For each one, estimate the corresponding
|
||||||
// scalar variable x as d-rS, with S the separator (remaining columns).
|
// scalar variable x as d-rS, with S the separator (remaining columns).
|
||||||
// Then update A and b by substituting x with d-rS, zero-ing out x's column.
|
// Then update A and b by substituting x with d-rS, zero-ing out x's column.
|
||||||
gtsam::print(Ab, "Ab = ");
|
if (verbose) gtsam::print(Ab, "Ab = ");
|
||||||
cout << " n = " << n << endl;
|
|
||||||
for (size_t j=0; j<n; ++j) {
|
for (size_t j=0; j<n; ++j) {
|
||||||
cout << "--------------------" << endl;
|
if (verbose) cout << "j: " << j << endl;
|
||||||
cout << "j: " << j << endl;
|
|
||||||
// extract the first column of A
|
// extract the first column of A
|
||||||
Vector a = Ab.col(j);
|
Vector a = Ab.col(j);
|
||||||
gtsam::print(a, "a = ");
|
|
||||||
|
|
||||||
// Calculate weighted pseudo-inverse and corresponding precision
|
// Calculate weighted pseudo-inverse and corresponding precision
|
||||||
gttic(constrained_QR_weightedPseudoinverse);
|
gttic(constrained_QR_weightedPseudoinverse);
|
||||||
double precision = weightedPseudoinverse(a, weights, pseudo);
|
double precision = weightedPseudoinverse(a, weights, pseudo);
|
||||||
gttoc(constrained_QR_weightedPseudoinverse);
|
gttoc(constrained_QR_weightedPseudoinverse);
|
||||||
cout << "precision: " << precision << endl;
|
if (verbose) gtsam::print(pseudo, "pseudo: ");
|
||||||
gtsam::print(pseudo, "pseudo: ");
|
|
||||||
|
|
||||||
// If precision is zero, no information on this column
|
// If precision is zero, no information on this column
|
||||||
// This is actually not limited to constraints, could happen in Gaussian::QR
|
// This is actually not limited to constraints, could happen in Gaussian::QR
|
||||||
|
|
@ -468,14 +463,12 @@ SharedDiagonal Constrained::QR(Matrix& Ab) const {
|
||||||
rd(j)=1.0; // put 1 on diagonal
|
rd(j)=1.0; // put 1 on diagonal
|
||||||
for (size_t j2=j+1; j2<n+1; ++j2) { // and fill in remainder with dot-products
|
for (size_t j2=j+1; j2<n+1; ++j2) { // and fill in remainder with dot-products
|
||||||
Vector Abj2 = Ab.col(j2);
|
Vector Abj2 = Ab.col(j2);
|
||||||
gtsam::print(Abj2, "Ab.col(j2): ");
|
|
||||||
rd(j2) = pseudo.dot(Ab.col(j2));
|
rd(j2) = pseudo.dot(Ab.col(j2));
|
||||||
}
|
}
|
||||||
gttoc(constrained_QR_create_rd);
|
gttoc(constrained_QR_create_rd);
|
||||||
|
|
||||||
// construct solution (r, d, sigma)
|
// construct solution (r, d, sigma)
|
||||||
Rd.push_back(boost::make_tuple(j, rd, precision));
|
Rd.push_back(boost::make_tuple(j, rd, precision));
|
||||||
gtsam::print(rd, "rd = ");
|
|
||||||
|
|
||||||
// exit after rank exhausted
|
// exit after rank exhausted
|
||||||
if (Rd.size()>=maxRank) break;
|
if (Rd.size()>=maxRank) break;
|
||||||
|
|
@ -484,16 +477,12 @@ SharedDiagonal Constrained::QR(Matrix& Ab) const {
|
||||||
gttic(constrained_QR_update_Ab);
|
gttic(constrained_QR_update_Ab);
|
||||||
Ab.middleCols(j+1,n-j) -= a * rd.segment(j+1, n-j).transpose();
|
Ab.middleCols(j+1,n-j) -= a * rd.segment(j+1, n-j).transpose();
|
||||||
gttoc(constrained_QR_update_Ab);
|
gttoc(constrained_QR_update_Ab);
|
||||||
gtsam::print(Ab, "Updated Ab = ");
|
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
cout << "Create storage for precisions" << endl;
|
|
||||||
// Create storage for precisions
|
// Create storage for precisions
|
||||||
Vector precisions(Rd.size());
|
Vector precisions(Rd.size());
|
||||||
|
|
||||||
gttic(constrained_QR_write_back_into_Ab);
|
gttic(constrained_QR_write_back_into_Ab);
|
||||||
cout << "write back result" << endl;
|
|
||||||
// Write back result in Ab, imperative as we are
|
// Write back result in Ab, imperative as we are
|
||||||
// TODO: test that is correct if a column was skipped !!!!
|
// TODO: test that is correct if a column was skipped !!!!
|
||||||
size_t i = 0; // start with first row
|
size_t i = 0; // start with first row
|
||||||
|
|
@ -511,8 +500,6 @@ SharedDiagonal Constrained::QR(Matrix& Ab) const {
|
||||||
}
|
}
|
||||||
gttoc(constrained_QR_write_back_into_Ab);
|
gttoc(constrained_QR_write_back_into_Ab);
|
||||||
|
|
||||||
cout << "return " << (int) mixed << endl;
|
|
||||||
gtsam::print(precisions, "precisions:");
|
|
||||||
// Must include mu, as the defaults might be higher, resulting in non-convergence
|
// Must include mu, as the defaults might be higher, resulting in non-convergence
|
||||||
return mixed ? Constrained::MixedPrecisions(mu_, precisions) : Diagonal::Precisions(precisions);
|
return mixed ? Constrained::MixedPrecisions(mu_, precisions) : Diagonal::Precisions(precisions);
|
||||||
}
|
}
|
||||||
|
|
|
||||||
|
|
@ -210,7 +210,7 @@ bool QPSolver::updateWorkingSetInplace(GaussianFactorGraph& workingGraph,
|
||||||
/* ************************************************************************* */
|
/* ************************************************************************* */
|
||||||
boost::tuple<double, int, int> QPSolver::computeStepSize(const GaussianFactorGraph& workingGraph,
|
boost::tuple<double, int, int> QPSolver::computeStepSize(const GaussianFactorGraph& workingGraph,
|
||||||
const VectorValues& xk, const VectorValues& p) const {
|
const VectorValues& xk, const VectorValues& p) const {
|
||||||
static bool debug = true;
|
static bool debug = false;
|
||||||
|
|
||||||
double minAlpha = 1.0;
|
double minAlpha = 1.0;
|
||||||
int closestFactorIx = -1, closestSigmaIx = -1;
|
int closestFactorIx = -1, closestSigmaIx = -1;
|
||||||
|
|
@ -228,9 +228,7 @@ boost::tuple<double, int, int> QPSolver::computeStepSize(const GaussianFactorGra
|
||||||
Vector aj = jacobian->getA(xj).row(s);
|
Vector aj = jacobian->getA(xj).row(s);
|
||||||
ajTp += aj.dot(pj);
|
ajTp += aj.dot(pj);
|
||||||
}
|
}
|
||||||
if (debug) {
|
if (debug) cout << "s, ajTp: " << s << " " << ajTp << endl;
|
||||||
cout << "s, ajTp: " << s << " " << ajTp << endl;
|
|
||||||
}
|
|
||||||
|
|
||||||
// Check if aj'*p >0. Don't care if it's not.
|
// Check if aj'*p >0. Don't care if it's not.
|
||||||
if (ajTp<=0) continue;
|
if (ajTp<=0) continue;
|
||||||
|
|
@ -242,15 +240,11 @@ boost::tuple<double, int, int> QPSolver::computeStepSize(const GaussianFactorGra
|
||||||
Vector aj = jacobian->getA(xj).row(s);
|
Vector aj = jacobian->getA(xj).row(s);
|
||||||
ajTx += aj.dot(xkj);
|
ajTx += aj.dot(xkj);
|
||||||
}
|
}
|
||||||
if (debug) {
|
if (debug) cout << "b[s], ajTx: " << b[s] << " " << ajTx << " " << ajTp << endl;
|
||||||
cout << "b[s], ajTx: " << b[s] << " " << ajTx << " " << ajTp << endl;
|
|
||||||
}
|
|
||||||
|
|
||||||
// alpha = (bj - aj'*xk) / (aj'*p)
|
// alpha = (bj - aj'*xk) / (aj'*p)
|
||||||
double alpha = (b[s] - ajTx)/ajTp;
|
double alpha = (b[s] - ajTx)/ajTp;
|
||||||
if (debug) {
|
if (debug) cout << "alpha: " << alpha << endl;
|
||||||
cout << "alpha: " << alpha << endl;
|
|
||||||
}
|
|
||||||
|
|
||||||
// We want the minimum of all those max alphas
|
// We want the minimum of all those max alphas
|
||||||
if (alpha < minAlpha) {
|
if (alpha < minAlpha) {
|
||||||
|
|
@ -266,7 +260,7 @@ boost::tuple<double, int, int> QPSolver::computeStepSize(const GaussianFactorGra
|
||||||
|
|
||||||
/* ************************************************************************* */
|
/* ************************************************************************* */
|
||||||
bool QPSolver::iterateInPlace(GaussianFactorGraph& workingGraph, VectorValues& currentSolution) const {
|
bool QPSolver::iterateInPlace(GaussianFactorGraph& workingGraph, VectorValues& currentSolution) const {
|
||||||
static bool debug = true;
|
static bool debug = false;
|
||||||
// Obtain the solution from the current working graph
|
// Obtain the solution from the current working graph
|
||||||
VectorValues newSolution = workingGraph.optimize();
|
VectorValues newSolution = workingGraph.optimize();
|
||||||
if (debug) newSolution.print("New solution:");
|
if (debug) newSolution.print("New solution:");
|
||||||
|
|
@ -294,9 +288,7 @@ bool QPSolver::iterateInPlace(GaussianFactorGraph& workingGraph, VectorValues& c
|
||||||
int factorIx, sigmaIx;
|
int factorIx, sigmaIx;
|
||||||
VectorValues p = newSolution - currentSolution;
|
VectorValues p = newSolution - currentSolution;
|
||||||
boost::tie(alpha, factorIx, sigmaIx) = computeStepSize(workingGraph, currentSolution, p);
|
boost::tie(alpha, factorIx, sigmaIx) = computeStepSize(workingGraph, currentSolution, p);
|
||||||
if (debug) {
|
if (debug) cout << "alpha, factorIx, sigmaIx: " << alpha << " " << factorIx << " " << sigmaIx << endl;
|
||||||
cout << "alpha, factorIx, sigmaIx: " << alpha << " " << factorIx << " " << sigmaIx << endl;
|
|
||||||
}
|
|
||||||
// also add to the working set the one that complains the most
|
// also add to the working set the one that complains the most
|
||||||
updateWorkingSetInplace(workingGraph, factorIx, sigmaIx, 0.0);
|
updateWorkingSetInplace(workingGraph, factorIx, sigmaIx, 0.0);
|
||||||
// step!
|
// step!
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue