From 2b2a20d9d90be356b3afdfbd5884578a686e5a16 Mon Sep 17 00:00:00 2001 From: mosfet80 Date: Mon, 6 Nov 2023 21:35:59 +0100 Subject: [PATCH] [planegcs] Removed unused code. (#10684) * Revert cleanplanegcs: removed unused code, removed redefinition of pi * Sketcher: Switch pi refs to double and constexpr * Modify code to use the new pi constant immediately --------- Co-authored-by: Chris Hennes --- src/Mod/Sketcher/App/planegcs/Constraints.cpp | 110 ++--------------- src/Mod/Sketcher/App/planegcs/GCS.cpp | 116 +----------------- src/Mod/Sketcher/App/planegcs/Geo.cpp | 21 ---- src/Mod/Sketcher/App/planegcs/Geo.h | 7 +- src/Mod/Sketcher/App/planegcs/SubSystem.cpp | 17 --- src/Mod/Sketcher/App/planegcs/Util.h | 5 - 6 files changed, 16 insertions(+), 260 deletions(-) diff --git a/src/Mod/Sketcher/App/planegcs/Constraints.cpp b/src/Mod/Sketcher/App/planegcs/Constraints.cpp index 7086ab4c28..aee8375afa 100644 --- a/src/Mod/Sketcher/App/planegcs/Constraints.cpp +++ b/src/Mod/Sketcher/App/planegcs/Constraints.cpp @@ -752,7 +752,6 @@ double ConstraintP2PDistance::grad(double* param) double ConstraintP2PDistance::maxStep(MAP_pD_D& dir, double lim) { MAP_pD_D::iterator it; - // distance() >= 0 it = dir.find(distance()); if (it != dir.end()) { if (it->second < 0.) { @@ -863,12 +862,11 @@ double ConstraintP2PAngle::grad(double* param) double ConstraintP2PAngle::maxStep(MAP_pD_D& dir, double lim) { - // step(angle()) <= pi/18 = 10° MAP_pD_D::iterator it = dir.find(angle()); if (it != dir.end()) { double step = std::abs(it->second); - if (step > M_PI / 18.0) { - lim = std::min(lim, (M_PI / 18.0) / step); + if (step > pi_18) { + lim = std::min(lim, pi_18 / step); } } return lim; @@ -908,18 +906,14 @@ double ConstraintP2LDistance::error() double dx = x2 - x1; double dy = y2 - y1; double d = sqrt(dx * dx + dy * dy); // line length - double area = - std::abs(-x0 * dy + y0 * dx + x1 * y2 - - x2 * y1); // = x1y2 - x2y1 - x0y2 + x2y0 + x0y1 - x1y0 = 2*(triangle area) + double area = std::abs(-x0 * dy + y0 * dx + x1 * y2 - x2 * y1); return scale * (area / d - dist); } double ConstraintP2LDistance::grad(double* param) { double deriv = 0.; - // darea/dx0 = (y1-y2) darea/dy0 = (x2-x1) - // darea/dx1 = (y2-y0) darea/dy1 = (x0-x2) - // darea/dx2 = (y0-y1) darea/dy2 = (x1-x0) + if (param == p0x() || param == p0y() || param == p1x() || param == p1y() || param == p2x() || param == p2y()) { double x0 = *p0x(), x1 = *p1x(), x2 = *p2x(); @@ -961,7 +955,6 @@ double ConstraintP2LDistance::grad(double* param) double ConstraintP2LDistance::maxStep(MAP_pD_D& dir, double lim) { MAP_pD_D::iterator it; - // distance() >= 0 it = dir.find(distance()); if (it != dir.end()) { if (it->second < 0.) { @@ -1056,17 +1049,13 @@ double ConstraintPointOnLine::error() double dx = x2 - x1; double dy = y2 - y1; double d = sqrt(dx * dx + dy * dy); - double area = -x0 * dy + y0 * dx + x1 * y2 - - x2 * y1; // = x1y2 - x2y1 - x0y2 + x2y0 + x0y1 - x1y0 = 2*(triangle area) + double area = -x0 * dy + y0 * dx + x1 * y2 - x2 * y1; return scale * area / d; } double ConstraintPointOnLine::grad(double* param) { double deriv = 0.; - // darea/dx0 = (y1-y2) darea/dy0 = (x2-x1) - // darea/dx1 = (y2-y0) darea/dy1 = (x0-x2) - // darea/dx2 = (y0-y1) darea/dy2 = (x1-x0) if (param == p0x() || param == p0y() || param == p1x() || param == p1y() || param == p2x() || param == p2y()) { double x0 = *p0x(), x1 = *p1x(), x2 = *p2x(); @@ -1453,12 +1442,11 @@ double ConstraintL2LAngle::grad(double* param) double ConstraintL2LAngle::maxStep(MAP_pD_D& dir, double lim) { - // step(angle()) <= pi/18 = 10° MAP_pD_D::iterator it = dir.find(angle()); if (it != dir.end()) { double step = std::abs(it->second); - if (step > M_PI / 18.0) { - lim = std::min(lim, (M_PI / 18.0) / step); + if (step > pi_18) { + lim = std::min(lim, pi_18 / step); } } return lim; @@ -1517,17 +1505,13 @@ double ConstraintMidpointOnLine::error() double dx = x2 - x1; double dy = y2 - y1; double d = sqrt(dx * dx + dy * dy); - double area = -x0 * dy + y0 * dx + x1 * y2 - - x2 * y1; // = x1y2 - x2y1 - x0y2 + x2y0 + x0y1 - x1y0 = 2*(triangle area) + double area = -x0 * dy + y0 * dx + x1 * y2 - x2 * y1; // = 2*(triangle area) return scale * area / d; } double ConstraintMidpointOnLine::grad(double* param) { double deriv = 0.; - // darea/dx0 = (y1-y2) darea/dy0 = (x2-x1) - // darea/dx1 = (y2-y0) darea/dy1 = (x0-x2) - // darea/dx2 = (y0-y1) darea/dy2 = (x1-x0) if (param == l1p1x() || param == l1p1y() || param == l1p2x() || param == l1p2y() || param == l2p1x() || param == l2p1y() || param == l2p2x() || param == l2p2y()) { double x0 = ((*l1p1x()) + (*l1p2x())) / 2; @@ -1831,23 +1815,6 @@ double ConstraintEllipseTangentLine::grad(double* param) double deriv; errorgrad(nullptr, &deriv, param); -// use numeric for testing -#if 0 - double const eps = 0.00001; - double oldparam = *param; - double v0 = this->error(); - *param += eps; - double vr = this->error(); - *param = oldparam - eps; - double vl = this->error(); - *param = oldparam; - //If not nasty, real derivative should be between left one and right one - double numretl = (v0 - vl) / eps; - double numretr = (vr - v0) / eps; - assert(deriv <= std::max(numretl, numretr)); - assert(deriv >= std::min(numretl, numretr)); -#endif - return deriv * scale; } @@ -1970,23 +1937,6 @@ double ConstraintInternalAlignmentPoint2Ellipse::grad(double* param) double deriv; errorgrad(nullptr, &deriv, param); -// use numeric for testing -#if 0 - double const eps = 0.00001; - double oldparam = *param; - double v0 = this->error(); - *param += eps; - double vr = this->error(); - *param = oldparam - eps; - double vl = this->error(); - *param = oldparam; - //If not nasty, real derivative should be between left one and right one - double numretl = (v0 - vl) / eps; - double numretr = (vr - v0) / eps; - assert(deriv <= std::max(numretl, numretr)); - assert(deriv >= std::min(numretl, numretr)); -#endif - return deriv * scale; } @@ -2374,15 +2324,6 @@ double ConstraintCurveValue::grad(double* param) double ConstraintCurveValue::maxStep(MAP_pD_D& /*dir*/, double lim) { - // step(angle()) <= pi/18 = 10° - /* TODO: curve-dependent parameter change limiting?? - MAP_pD_D::iterator it = dir.find(this->u()); - if (it != dir.end()) { - double step = std::abs(it->second); - if (step > M_PI/18.) - lim = std::min(lim, (M_PI/18.) / step); - } - */ return lim; } @@ -2706,24 +2647,6 @@ double ConstraintAngleViaPoint::grad(double* param) deriv -= ((-n1.dx) * n1.y / pow(n1.length(), 2) + n1.dy * n1.x / pow(n1.length(), 2)); deriv += ((-n2.dx) * n2.y / pow(n2.length(), 2) + n2.dy * n2.x / pow(n2.length(), 2)); - -// use numeric for testing -#if 0 - double const eps = 0.00001; - double oldparam = *param; - double v0 = this->error(); - *param += eps; - double vr = this->error(); - *param = oldparam - eps; - double vl = this->error(); - *param = oldparam; - //If not nasty, real derivative should be between left one and right one - double numretl = (v0-vl)/eps; - double numretr = (vr-v0)/eps; - assert(deriv <= std::max(numretl,numretr) ); - assert(deriv >= std::min(numretl,numretr) ); -#endif - return scale * deriv; } @@ -2841,23 +2764,6 @@ double ConstraintSnell::grad(double* param) double deriv; errorgrad(nullptr, &deriv, param); -// use numeric for testing -#if 0 - double const eps = 0.00001; - double oldparam = *param; - double v0 = this->error(); - *param += eps; - double vr = this->error(); - *param = oldparam - eps; - double vl = this->error(); - *param = oldparam; - //If not nasty, real derivative should be between left one and right one - double numretl = (v0 - vl) / eps; - double numretr = (vr - v0) / eps; - assert(deriv <= std::max(numretl, numretr)); - assert(deriv >= std::min(numretl, numretr)); -#endif - return scale * deriv; } diff --git a/src/Mod/Sketcher/App/planegcs/GCS.cpp b/src/Mod/Sketcher/App/planegcs/GCS.cpp index 19d59ecae5..1b7bf211bc 100644 --- a/src/Mod/Sketcher/App/planegcs/GCS.cpp +++ b/src/Mod/Sketcher/App/planegcs/GCS.cpp @@ -26,10 +26,6 @@ #pragma warning(disable : 4996) #endif -// #define _GCS_DEBUG -// #define _GCS_DEBUG_SOLVER_JACOBIAN_QR_DECOMPOSITION_TRIANGULAR_MATRIX -// #define _DEBUG_TO_FILE // Many matrices surpass the report view string size. -// #define PROFILE_DIAGNOSE #undef _GCS_DEBUG #undef _GCS_DEBUG_SOLVER_JACOBIAN_QR_DECOMPOSITION_TRIANGULAR_MATRIX #undef _DEBUG_TO_FILE @@ -505,76 +501,6 @@ System::System() #endif } -/*DeepSOIC: seriously outdated, needs redesign -System::System(std::vector clist_) -: plist(0), - c2p(), p2c(), - subSystems(0), subSystemsAux(0), - reference(0), - hasUnknowns(false), hasDiagnosis(false), isInit(false) -{ - // create own (shallow) copy of constraints - for (std::vector::iterator constr=clist_.begin(); - constr != clist_.end(); ++constr) { - Constraint *newconstr = 0; - switch ((*constr)->getTypeId()) { - case Equal: { - ConstraintEqual *oldconstr = static_cast(*constr); - newconstr = new ConstraintEqual(*oldconstr); - break; - } - case Difference: { - ConstraintDifference *oldconstr = static_cast(*constr); - newconstr = new ConstraintDifference(*oldconstr); - break; - } - case P2PDistance: { - ConstraintP2PDistance *oldconstr = static_cast(*constr); - newconstr = new ConstraintP2PDistance(*oldconstr); - break; - } - case P2PAngle: { - ConstraintP2PAngle *oldconstr = static_cast(*constr); - newconstr = new ConstraintP2PAngle(*oldconstr); - break; - } - case P2LDistance: { - ConstraintP2LDistance *oldconstr = static_cast(*constr); - newconstr = new ConstraintP2LDistance(*oldconstr); - break; - } - case PointOnLine: { - ConstraintPointOnLine *oldconstr = static_cast(*constr); - newconstr = new ConstraintPointOnLine(*oldconstr); - break; - } - case Parallel: { - ConstraintParallel *oldconstr = static_cast(*constr); - newconstr = new ConstraintParallel(*oldconstr); - break; - } - case Perpendicular: { - ConstraintPerpendicular *oldconstr = static_cast(*constr); newconstr = new ConstraintPerpendicular(*oldconstr); break; - } - case L2LAngle: { - ConstraintL2LAngle *oldconstr = static_cast(*constr); - newconstr = new ConstraintL2LAngle(*oldconstr); - break; - } - case MidpointOnLine: { - ConstraintMidpointOnLine *oldconstr = static_cast(*constr); newconstr = new ConstraintMidpointOnLine(*oldconstr); break; - } - case None: - break; - } - if (newconstr) - addConstraint(newconstr); - } -} -*/ - System::~System() { clear(); @@ -1057,7 +983,7 @@ int System::addConstraintPerpendicularLine2Arc(Point& p1, return addConstraintP2PAngle(p1, p2, a.startAngle, 0, tagId, driving); } else { - return addConstraintP2PAngle(p1, p2, a.startAngle, M_PI, tagId, driving); + return addConstraintP2PAngle(p1, p2, a.startAngle, pi, tagId, driving); } } @@ -1074,7 +1000,7 @@ int System::addConstraintPerpendicularArc2Line(Arc& a, return addConstraintP2PAngle(p1, p2, a.endAngle, 0, tagId, driving); } else { - return addConstraintP2PAngle(p1, p2, a.endAngle, M_PI, tagId, driving); + return addConstraintP2PAngle(p1, p2, a.endAngle, pi, tagId, driving); } } @@ -1085,7 +1011,7 @@ int System::addConstraintPerpendicularCircle2Arc(Point& center, bool driving) { addConstraintP2PDistance(a.start, center, radius, tagId, driving); - double incrAngle = *(a.startAngle) < *(a.endAngle) ? M_PI / 2 : -M_PI / 2; + double incrAngle = *(a.startAngle) < *(a.endAngle) ? pi_2 : -pi_2; double tangAngle = *a.startAngle + incrAngle; double dx = *(a.start.x) - *(center.x); double dy = *(a.start.y) - *(center.y); @@ -1104,7 +1030,7 @@ int System::addConstraintPerpendicularArc2Circle(Arc& a, bool driving) { addConstraintP2PDistance(a.end, center, radius, tagId, driving); - double incrAngle = *(a.startAngle) < *(a.endAngle) ? -M_PI / 2 : M_PI / 2; + double incrAngle = *(a.startAngle) < *(a.endAngle) ? -pi_2 : pi_2; double tangAngle = *a.endAngle + incrAngle; double dx = *(a.end.x) - *(center.x); double dy = *(a.end.y) - *(center.y); @@ -1331,14 +1257,6 @@ int System::addConstraintInternalAlignmentEllipseMajorDiameter(Ellipse& e, double Y_F1 = *e.focus1.y; double b = *e.radmin; - // P1=vector([X_1,Y_1]) - // P2=vector([X_2,Y_2]) - // dF1= (F1-C)/sqrt((F1-C)*(F1-C)) - // print "these are the extreme points of the major axis" - // PA = C + a * dF1 - // PN = C - a * dF1 - // print "this is a simple function to know which point is closer to the positive edge of the - // ellipse" DMC=(P1-PA)*(P1-PA)-(P2-PA)*(P2-PA) double closertopositivemajor = pow(X_1 - X_c - (X_F1 - X_c) * sqrt(pow(b, 2) + pow(X_F1 - X_c, 2) + pow(Y_F1 - Y_c, 2)) @@ -1397,8 +1315,6 @@ int System::addConstraintInternalAlignmentEllipseMinorDiameter(Ellipse& e, double Y_F1 = *e.focus1.y; double b = *e.radmin; - // Same idea as for major above, but for minor - // DMC=(P1-PA)*(P1-PA)-(P2-PA)*(P2-PA) double closertopositiveminor = pow(X_1 - X_c + b * (Y_F1 - Y_c) / sqrt(pow(X_F1 - X_c, 2) + pow(Y_F1 - Y_c, 2)), 2) - pow(X_2 - X_c + b * (Y_F1 - Y_c) / sqrt(pow(X_F1 - X_c, 2) + pow(Y_F1 - Y_c, 2)), 2) @@ -1465,14 +1381,6 @@ int System::addConstraintInternalAlignmentHyperbolaMajorDiameter(Hyperbola& e, double Y_F1 = *e.focus1.y; double b = *e.radmin; - // P1=vector([X_1,Y_1]) - // P2=vector([X_2,Y_2]) - // dF1= (F1-C)/sqrt((F1-C)*(F1-C)) - // print "these are the extreme points of the major axis" - // PA = C + a * dF1 - // PN = C - a * dF1 - // print "this is a simple function to know which point is closer to the positive edge of the - // ellipse" DMC=(P1-PA)*(P1-PA)-(P2-PA)*(P2-PA) double closertopositivemajor = pow(-X_1 + X_c + (X_F1 - X_c) * (-pow(b, 2) + pow(X_F1 - X_c, 2) + pow(Y_F1 - Y_c, 2)) @@ -1555,8 +1463,6 @@ int System::addConstraintInternalAlignmentHyperbolaMinorDiameter(Hyperbola& e, double Y_F1 = *e.focus1.y; double b = *e.radmin; - // Same idea as for major above, but for minor - // DMC=(P1-PA)*(P1-PA)-(P2-PA)*(P2-PA) double closertopositiveminor = pow(-X_1 + X_c + b * (Y_F1 - Y_c) / sqrt(pow(X_F1 - X_c, 2) + pow(Y_F1 - Y_c, 2)) + (X_F1 - X_c) * (-pow(b, 2) + pow(X_F1 - X_c, 2) + pow(Y_F1 - Y_c, 2)) @@ -2505,11 +2411,6 @@ int System::solve_DL(SubSystem* subsys, bool isRedundantsolving) break; } - // it didn't work in some tests - // // restrict h_dl according to maxStep - // double scale = subsys->maxStep(h_dl); - // if (scale < 1.) - // h_dl *= scale; // get the new values double err_new; @@ -4794,12 +4695,6 @@ int System::solve(SubSystem* subsysA, SubSystem* subsysB, bool /*isFine*/, bool double alpha = 1; alpha = std::min(alpha, subsysA->maxStep(plistAB, xdir)); - // From the book "Numerical Optimization - Jorge Nocedal, Stephen J. Wright". - // See https://forum.freecad.org/viewtopic.php?f=10&t=35469. - // Eq. 18.32 - // double mu = lambda.lpNorm() + 0.01; - // Eq. 18.33 - // double mu = grad.dot(xdir) / ( (1.-rho) * resA.lpNorm<1>()); // Eq. 18.36 mu = std::max(mu, (grad.dot(xdir) + std::max(0., 0.5 * xdir.dot(B * xdir))) @@ -4821,9 +4716,6 @@ int System::solve(SubSystem* subsysA, SubSystem* subsysB, bool /*isFine*/, bool bool first = true; while (f > f0 + eta * alpha * deriv) { if (first) { - // try a second order step - // xdir1 = JA.jacobiSvd(Eigen::ComputeThinU | - // Eigen::ComputeThinV).solve(-resA); xdir1 = -Y * resA; x += xdir1; // = x0 + alpha * xdir + xdir1 subsysA->setParams(plistAB, x); diff --git a/src/Mod/Sketcher/App/planegcs/Geo.cpp b/src/Mod/Sketcher/App/planegcs/Geo.cpp index cba43d830a..b085384fd6 100644 --- a/src/Mod/Sketcher/App/planegcs/Geo.cpp +++ b/src/Mod/Sketcher/App/planegcs/Geo.cpp @@ -325,27 +325,6 @@ DeriVector2 Ellipse::CalculateNormal(const Point& p, const double* derivparam) c // return sum of normalized pf2, pf2 DeriVector2 ret = pf1.getNormalized().sum(pf2.getNormalized()); -// numeric derivatives for testing -#if 0 // make sure to enable DEBUG_DERIVS when enabling - if(derivparam) { - double const eps = 0.00001; - double oldparam = *derivparam; - DeriVector2 v0 = this->CalculateNormal(p); - *derivparam += eps; - DeriVector2 vr = this->CalculateNormal(p); - *derivparam = oldparam - eps; - DeriVector2 vl = this->CalculateNormal(p); - *derivparam = oldparam; - //If not nasty, real derivative should be between left one and right one - DeriVector2 numretl ((v0.x-vl.x)/eps, (v0.y-vl.y)/eps); - DeriVector2 numretr ((vr.x-v0.x)/eps, (vr.y-v0.y)/eps); - assert(ret.dx <= std::max(numretl.x,numretr.x) ); - assert(ret.dx >= std::min(numretl.x,numretr.x) ); - assert(ret.dy <= std::max(numretl.y,numretr.y) ); - assert(ret.dy >= std::min(numretl.y,numretr.y) ); - } -#endif - return ret; } diff --git a/src/Mod/Sketcher/App/planegcs/Geo.h b/src/Mod/Sketcher/App/planegcs/Geo.h index a6f6487274..91f0eb6f6f 100644 --- a/src/Mod/Sketcher/App/planegcs/Geo.h +++ b/src/Mod/Sketcher/App/planegcs/Geo.h @@ -23,10 +23,8 @@ #ifndef PLANEGCS_GEO_H #define PLANEGCS_GEO_H -#include - #include "Util.h" - +#include namespace GCS { @@ -50,6 +48,9 @@ public: }; using VEC_P = std::vector; +static constexpr double pi = boost::math::constants::pi(); +static constexpr double pi_2 = pi / 2.0; +static constexpr double pi_18 = pi / 18.0; /// Class DeriVector2 holds a vector value and its derivative on the /// parameter that the derivatives are being calculated for now. x,y is the diff --git a/src/Mod/Sketcher/App/planegcs/SubSystem.cpp b/src/Mod/Sketcher/App/planegcs/SubSystem.cpp index 494d6a8fd4..9d88131681 100644 --- a/src/Mod/Sketcher/App/planegcs/SubSystem.cpp +++ b/src/Mod/Sketcher/App/planegcs/SubSystem.cpp @@ -258,22 +258,6 @@ void SubSystem::calcResidual(Eigen::VectorXd& r, double& err) err *= 0.5; } -/* -void SubSystem::calcJacobi() -{ - assert(grad.size() != xsize); - - for (MAP_pD_pD::const_iterator param=pmap.begin(); - param != pmap.end(); ++param) { - // assert(p2c.find(param->second) != p2c.end()); - std::vector constrs=p2c[param->second]; - for (std::vector::const_iterator constr = constrs.begin(); - constr != constrs.end(); ++constr) - jacobi.set(*constr,param->second,(*constr)->grad(param->second)); - } -} -*/ - void SubSystem::calcJacobi(VEC_pD& params, Eigen::MatrixXd& jacobi) { jacobi.setZero(csize, params.size()); @@ -300,7 +284,6 @@ void SubSystem::calcGrad(VEC_pD& params, Eigen::VectorXd& grad) for (int j = 0; j < int(params.size()); j++) { MAP_pD_pD::const_iterator pmapfind = pmap.find(params[j]); if (pmapfind != pmap.end()) { - // assert(p2c.find(pmapfind->second) != p2c.end()); std::vector constrs = p2c[pmapfind->second]; for (std::vector::const_iterator constr = constrs.begin(); constr != constrs.end(); diff --git a/src/Mod/Sketcher/App/planegcs/Util.h b/src/Mod/Sketcher/App/planegcs/Util.h index 644b5a7e3c..87b0892775 100644 --- a/src/Mod/Sketcher/App/planegcs/Util.h +++ b/src/Mod/Sketcher/App/planegcs/Util.h @@ -38,11 +38,6 @@ using MAP_pD_D = std::map; using MAP_pD_I = std::map; using SET_pD = std::set; using SET_I = std::set; - -#ifndef M_PI -#define M_PI 3.14159265358979323846 -#endif - } // namespace GCS #endif // PLANEGCS_UTIL_H