diff --git a/src/system.cpp b/src/system.cpp index 89d315938..ded9e9eab 100644 --- a/src/system.cpp +++ b/src/system.cpp @@ -291,13 +291,59 @@ bool System::SolveLeastSquares() { } } - SparseMatrix AAt = mat.A.num * mat.A.num.transpose(); - AAt.makeCompressed(); - VectorXd z(mat.n); + if(mat.m == 0 || mat.n == 0) { + mat.X = VectorXd::Zero(mat.n); + return true; + } + + if(mat.m <= mat.n) { + // We want the minimum norm solution of the underdetermined system + // A*x = B, which is x = A'*(A*A')^-1*B. Forming the normal equations + // A*A' explicitly squares the condition number of A; since our + // Jacobian mixes dimensionless quantities with lengths, areas and + // volumes, its condition number already grows with the size of the + // sketch, and squaring it is what makes large models fail. Eigen's + // rank-revealing QR then drops pivots below a threshold proportional + // to the largest column norm, and its rank-truncated solve silently + // zeroes the corresponding component of the Newton step, so the + // residual of one equation can never be driven to zero and we report + // unsolvable constraints for a perfectly solvable sketch. + // + // Factor A' instead, as A'*P = Q*R. Then A*A' = P*R'*R*P', so + // z = (A*A')^-1*B = P*R^-1*R'^-1*P'*B, + // which is two triangular solves, and the rank decision is taken on + // R (whose pivots scale like A, not like A*A'). + SparseMatrix At = mat.A.num.transpose(); + At.makeCompressed(); + + SparseQR, COLAMDOrdering> solver; + solver.compute(At); + if(solver.info() != Success) return false; + + const int rank = (int)solver.rank(); + VectorXd rhs = solver.colsPermutation().transpose() * mat.B.num; + VectorXd v = VectorXd::Zero(mat.m); + if(rank > 0) { + SparseMatrix R = solver.matrixR().topLeftCorner(rank, rank); + VectorXd w = R.triangularView().transpose() + .solve(rhs.topRows(rank)); + v.topRows(rank) = R.triangularView().solve(w); + } + VectorXd z = solver.colsPermutation() * v; - if(!SolveLinearSystem(AAt, mat.B.num, &z)) return false; + mat.X = mat.A.num.transpose() * z; + } else { + // More equations than unknowns; A' is wide, which Eigen's sparse QR + // does not handle, so fall back to the normal equations. A system + // like that is redundant anyway, and gets reported as such. + SparseMatrix AAt = mat.A.num * mat.A.num.transpose(); + AAt.makeCompressed(); + VectorXd z(mat.m); + + if(!SolveLinearSystem(AAt, mat.B.num, &z)) return false; - mat.X = mat.A.num.transpose() * z; + mat.X = mat.A.num.transpose() * z; + } for(int c = 0; c < mat.n; c++) { mat.X[c] *= scale[c]; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 0dafd42db..51d1cbd80 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -54,6 +54,7 @@ set(testsuite_SOURCES constraint/equal_radius/test.cpp constraint/where_dragged/test.cpp constraint/comment/test.cpp + constraint/large_dimensions/test.cpp request/arc_of_circle/test.cpp request/circle/test.cpp request/cubic/test.cpp diff --git a/test/constraint/large_dimensions/perpendicular_4m.slvs b/test/constraint/large_dimensions/perpendicular_4m.slvs new file mode 100644 index 000000000..bde75e04f --- /dev/null +++ b/test/constraint/large_dimensions/perpendicular_4m.slvs @@ -0,0 +1,381 @@ +±²³SolveSpaceREVa + + +Group.h.v=00000001 +Group.type=5000 +Group.name=#references +Group.color=ff000000 +Group.skipFirst=0 +Group.predef.swapUV=0 +Group.predef.negateU=0 +Group.predef.negateV=0 +Group.visible=1 +Group.suppress=0 +Group.relaxConstraints=0 +Group.allowRedundant=0 +Group.allDimsReference=0 +Group.remap={ +} +AddGroup + +Group.h.v=00000002 +Group.type=5001 +Group.order=1 +Group.name=sketch-in-plane +Group.activeWorkplane.v=80020000 +Group.color=ff000000 +Group.subtype=6000 +Group.skipFirst=0 +Group.predef.q.w=1.00000000000000000000 +Group.predef.origin.v=00010001 +Group.predef.swapUV=0 +Group.predef.negateU=0 +Group.predef.negateV=0 +Group.visible=1 +Group.suppress=0 +Group.relaxConstraints=0 +Group.allowRedundant=0 +Group.allDimsReference=0 +Group.remap={ +} +AddGroup + +Param.h.v.=00010010 +AddParam + +Param.h.v.=00010011 +AddParam + +Param.h.v.=00010012 +AddParam + +Param.h.v.=00010020 +Param.val=1.00000000000000000000 +AddParam + +Param.h.v.=00010021 +AddParam + +Param.h.v.=00010022 +AddParam + +Param.h.v.=00010023 +AddParam + +Param.h.v.=00020010 +AddParam + +Param.h.v.=00020011 +AddParam + +Param.h.v.=00020012 +AddParam + +Param.h.v.=00020020 +Param.val=0.50000000000000000000 +AddParam + +Param.h.v.=00020021 +Param.val=0.50000000000000000000 +AddParam + +Param.h.v.=00020022 +Param.val=0.50000000000000000000 +AddParam + +Param.h.v.=00020023 +Param.val=0.50000000000000000000 +AddParam + +Param.h.v.=00030010 +AddParam + +Param.h.v.=00030011 +AddParam + +Param.h.v.=00030012 +AddParam + +Param.h.v.=00030020 +Param.val=0.50000000000000000000 +AddParam + +Param.h.v.=00030021 +Param.val=-0.50000000000000000000 +AddParam + +Param.h.v.=00030022 +Param.val=-0.50000000000000000000 +AddParam + +Param.h.v.=00030023 +Param.val=-0.50000000000000000000 +AddParam + +Param.h.v.=00040010 +Param.val=-873.94336616523173688620 +AddParam + +Param.h.v.=00040011 +Param.val=888.35750558391430331540 +AddParam + +Param.h.v.=00040013 +Param.val=-310.42323046797423558019 +AddParam + +Param.h.v.=00040014 +Param.val=-3071.74919205236892594257 +AddParam + +Param.h.v.=00050010 +Param.val=-943.11295456966990968795 +AddParam + +Param.h.v.=00050011 +Param.val=-1154.79416602316041462473 +AddParam + +Param.h.v.=00050013 +Param.val=-4920.05203201180302130524 +AddParam + +Param.h.v.=00050014 +Param.val=-1583.69459073004077254154 +AddParam + +Param.h.v.=40000002 +Param.val=0.56702743975097524842 +AddParam + +Request.h.v=00000001 +Request.type=100 +Request.group.v=00000001 +Request.construction=0 +AddRequest + +Request.h.v=00000002 +Request.type=100 +Request.group.v=00000001 +Request.construction=0 +AddRequest + +Request.h.v=00000003 +Request.type=100 +Request.group.v=00000001 +Request.construction=0 +AddRequest + +Request.h.v=00000004 +Request.type=200 +Request.workplane.v=80020000 +Request.group.v=00000002 +Request.construction=0 +AddRequest + +Request.h.v=00000005 +Request.type=200 +Request.workplane.v=80020000 +Request.group.v=00000002 +Request.construction=0 +AddRequest + +Entity.h.v=00010000 +Entity.type=10000 +Entity.construction=0 +Entity.point[0].v=00010001 +Entity.normal.v=00010020 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00010001 +Entity.type=2000 +Entity.construction=1 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00010020 +Entity.type=3000 +Entity.construction=0 +Entity.point[0].v=00010001 +Entity.actNormal.w=1.00000000000000000000 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00020000 +Entity.type=10000 +Entity.construction=0 +Entity.point[0].v=00020001 +Entity.normal.v=00020020 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00020001 +Entity.type=2000 +Entity.construction=1 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00020020 +Entity.type=3000 +Entity.construction=0 +Entity.point[0].v=00020001 +Entity.actNormal.w=0.50000000000000000000 +Entity.actNormal.vx=0.50000000000000000000 +Entity.actNormal.vy=0.50000000000000000000 +Entity.actNormal.vz=0.50000000000000000000 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00030000 +Entity.type=10000 +Entity.construction=0 +Entity.point[0].v=00030001 +Entity.normal.v=00030020 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00030001 +Entity.type=2000 +Entity.construction=1 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00030020 +Entity.type=3000 +Entity.construction=0 +Entity.point[0].v=00030001 +Entity.actNormal.w=0.50000000000000000000 +Entity.actNormal.vx=-0.50000000000000000000 +Entity.actNormal.vy=-0.50000000000000000000 +Entity.actNormal.vz=-0.50000000000000000000 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00040000 +Entity.type=11000 +Entity.construction=0 +Entity.point[0].v=00040001 +Entity.point[1].v=00040002 +Entity.workplane.v=80020000 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00040001 +Entity.type=2001 +Entity.construction=0 +Entity.workplane.v=80020000 +Entity.actPoint.x=-873.94336616523173688620 +Entity.actPoint.y=888.35750558391430331540 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00040002 +Entity.type=2001 +Entity.construction=0 +Entity.workplane.v=80020000 +Entity.actPoint.x=-310.42323046797423558019 +Entity.actPoint.y=-3071.74919205236892594257 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00050000 +Entity.type=11000 +Entity.construction=0 +Entity.point[0].v=00050001 +Entity.point[1].v=00050002 +Entity.workplane.v=80020000 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00050001 +Entity.type=2001 +Entity.construction=0 +Entity.workplane.v=80020000 +Entity.actPoint.x=-943.11295456966990968795 +Entity.actPoint.y=-1154.79416602316041462473 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00050002 +Entity.type=2001 +Entity.construction=0 +Entity.workplane.v=80020000 +Entity.actPoint.x=-4920.05203201180302130524 +Entity.actPoint.y=-1583.69459073004077254154 +Entity.actVisible=1 +AddEntity + +Entity.h.v=80020000 +Entity.type=10000 +Entity.construction=0 +Entity.point[0].v=80020002 +Entity.normal.v=80020001 +Entity.actVisible=1 +AddEntity + +Entity.h.v=80020001 +Entity.type=3010 +Entity.construction=0 +Entity.point[0].v=80020002 +Entity.actNormal.w=1.00000000000000000000 +Entity.actVisible=1 +AddEntity + +Entity.h.v=80020002 +Entity.type=2012 +Entity.construction=1 +Entity.actVisible=1 +AddEntity + +Constraint.h.v=00000002 +Constraint.type=42 +Constraint.group.v=00000002 +Constraint.workplane.v=80020000 +Constraint.valP.v=40000002 +Constraint.ptA.v=00050001 +Constraint.entityA.v=00040000 +Constraint.other=0 +Constraint.other2=0 +Constraint.reference=0 +AddConstraint + +Constraint.h.v=00000003 +Constraint.type=30 +Constraint.group.v=00000002 +Constraint.workplane.v=80020000 +Constraint.valA=4000.00000000000000000000 +Constraint.ptA.v=00040001 +Constraint.ptB.v=00040002 +Constraint.other=0 +Constraint.other2=0 +Constraint.reference=0 +Constraint.disp.offset.x=497.85777038281713657852 +Constraint.disp.offset.y=-105.42002325418791031097 +AddConstraint + +Constraint.h.v=00000004 +Constraint.type=30 +Constraint.group.v=00000002 +Constraint.workplane.v=80020000 +Constraint.valA=4000.00000000000000000000 +Constraint.ptA.v=00050001 +Constraint.ptB.v=00050002 +Constraint.other=0 +Constraint.other2=0 +Constraint.reference=0 +Constraint.disp.offset.x=-27.53621872585551599855 +Constraint.disp.offset.y=264.80569179149119918293 +AddConstraint + +Constraint.h.v=00000005 +Constraint.type=122 +Constraint.group.v=00000002 +Constraint.workplane.v=80020000 +Constraint.entityA.v=00050000 +Constraint.entityB.v=00040000 +Constraint.other=0 +Constraint.other2=0 +Constraint.reference=0 +AddConstraint + diff --git a/test/constraint/large_dimensions/test.cpp b/test/constraint/large_dimensions/test.cpp new file mode 100644 index 000000000..a242d728e --- /dev/null +++ b/test/constraint/large_dimensions/test.cpp @@ -0,0 +1,42 @@ +#include "solvespace.h" + +#include "harness.h" + +// The two line segments of the sketch, each dimensioned to 4 m, constrained +// perpendicular to each other with the start of the second one on the first. +static const hEntity LINE_A_A = { 0x00040001 }, LINE_A_B = { 0x00040002 }; +static const hEntity LINE_B_A = { 0x00050001 }, LINE_B_B = { 0x00050002 }; + +static Vector Dir(hEntity a, hEntity b) { + return SK.GetEntity(b)->PointGetNum().Minus(SK.GetEntity(a)->PointGetNum()); +} + +// A sketch whose dimensions are large enough that the entries of the Jacobian +// span several orders of magnitude. Solving the Newton step through the normal +// equations squared that spread, which pushed a perfectly good pivot under the +// threshold below which a column counts as linearly dependent; the step then +// came back with that component zeroed, and the solver reported these entirely +// compatible constraints as incompatible. See #1354. +TEST_CASE(perpendicular_4m) { + // Unlike most fixtures in this suite, perpendicular_4m.slvs is the + // reporter's file exactly as it was uploaded, and it must stay that way. + // Re-saving it from SolveSpace writes out the *solved* parameter values, + // so the sketch loads already at its solution, the first Newton step is + // the zero step, and the bug this test is here for cannot happen. Do not + // canonicalize this file. + CHECK_LOAD("perpendicular_4m.slvs"); + + for(Group &g : SK.group) { + CHECK_TRUE(g.IsSolvedOkay()); + } + + // And it solved it correctly, not just to something. + Vector da = Dir(LINE_A_A, LINE_A_B), db = Dir(LINE_B_A, LINE_B_B); + CHECK_EQ_EPS(da.Magnitude(), 4000.0); + CHECK_EQ_EPS(db.Magnitude(), 4000.0); + CHECK_EQ_EPS(da.WithMagnitude(1).Dot(db.WithMagnitude(1)), 0.0); + + // The start of the second line lies on the first one. + Vector onLine = Dir(LINE_A_A, LINE_B_A); + CHECK_EQ_EPS(onLine.Cross(da).Magnitude() / da.Magnitude(), 0.0); +} diff --git a/test/debugtool.cpp b/test/debugtool.cpp index aea6b8475..41d973a88 100644 --- a/test/debugtool.cpp +++ b/test/debugtool.cpp @@ -4,15 +4,82 @@ // Copyright 2017 whitequark //----------------------------------------------------------------------------- +#include "solvespace.h" #include "expr.h" #include "platform/platform.h" using namespace SolveSpace; +static const char *SolveResultName(SolveResult how) { + switch(how) { + case SolveResult::OKAY: return "OKAY"; + case SolveResult::DIDNT_CONVERGE: return "DIDNT_CONVERGE"; + case SolveResult::REDUNDANT_OKAY: return "REDUNDANT_OKAY"; + case SolveResult::REDUNDANT_DIDNT_CONVERGE: return "REDUNDANT_DIDNT_CONVERGE"; + case SolveResult::TOO_MANY_UNKNOWNS: return "TOO_MANY_UNKNOWNS"; + } + return "?"; +} + +// Load a file, regenerate it, and report how every group solved. Useful to +// reproduce solver failures without a GUI. +static int CmdSolve(const std::string &filename) { + SS.Init(); + SS.showToolbar = false; + SS.checkClosedContour = false; + + if(!SS.LoadFromFile(Platform::Path::From(filename))) { + fprintf(stderr, "cannot load: %s\n", filename.c_str()); + return 1; + } + SS.AfterNewFile(); + + int failed = 0; + for(Group &g : SK.group) { + bool ok = g.IsSolvedOkay(); + if(!ok) failed++; + fprintf(stderr, "group %08x %-24s %-24s dof=%d%s\n", g.h.v, + g.DescriptionString().c_str(), SolveResultName(g.solved.how), + g.solved.dof, ok ? "" : " FAILED"); + for(int i = 0; i < g.solved.remove.n; i++) { + Constraint *c = SK.constraint.FindByIdNoOops(g.solved.remove[i]); + if(!c) continue; + fprintf(stderr, " bad constraint %08x %s\n", c->h.v, + c->DescriptionString().c_str()); + } + } + fprintf(stderr, "%s: %d group(s) failed to solve\n", filename.c_str(), failed); + return failed == 0 ? 0 : 1; +} + +// Load a file, regenerate it, and write it back out. This is what happens when +// a part that other files link to is opened, edited and saved again. +static int CmdResave(const std::string &filename) { + SS.Init(); + SS.showToolbar = false; + SS.checkClosedContour = false; + + Platform::Path path = Platform::Path::From(filename); + if(!SS.LoadFromFile(path)) { + fprintf(stderr, "cannot load: %s\n", filename.c_str()); + return 1; + } + SS.AfterNewFile(); + if(!SS.SaveToFile(path)) { + fprintf(stderr, "cannot save: %s\n", filename.c_str()); + return 1; + } + return 0; +} + int main(int argc, char **argv) { std::vector args = Platform::InitCli(argc, argv); - if(args.size() == 3 && args[1] == "expr") { + if(args.size() == 3 && args[1] == "solve") { + return CmdSolve(args[2]); + } else if(args.size() == 3 && args[1] == "resave") { + return CmdResave(args[2]); + } else if(args.size() == 3 && args[1] == "expr") { std::string expr = args[2], err; Expr *e = Expr::Parse(expr.c_str(), &err); if(e == NULL) { @@ -28,6 +95,10 @@ int main(int argc, char **argv) { Commands: expr [expr] Evaluate an expression. + solve [file.slvs] + Load a file and report how each group solved. + resave [file.slvs] + Load a file, regenerate it, and save it back. )"); }