kindred/create
1
1
Fork 0
Code Issues 30 Pull Requests Actions Packages Projects 9 Releases 2 Wiki Activity

[planegcs] Removed unused code. (#10684)

Browse Source 
* 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 <chennes@pioneerlibrarysystem.org>
This commit is contained in:
mosfet80
2023-11-06 21:35:59 +01:00
committed by GitHub
parent a7e1cba6b0
commit 2b2a20d9d9
6 changed files with 16 additions and 260 deletions
Download Patch File Download Diff File Expand all files Collapse all files

110
src/Mod/Sketcher/App/planegcs/Constraints.cpp
View File

@@ -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;
}

116
src/Mod/Sketcher/App/planegcs/GCS.cpp
View File

@@ -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<Constraint *> 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<Constraint *>::iterator constr=clist_.begin();
constr != clist_.end(); ++constr) {
Constraint *newconstr = 0;
switch ((*constr)->getTypeId()) {
case Equal: {
ConstraintEqual *oldconstr = static_cast<ConstraintEqual *>(*constr);
newconstr = new ConstraintEqual(*oldconstr);
break;
}
case Difference: {
ConstraintDifference *oldconstr = static_cast<ConstraintDifference *>(*constr);
newconstr = new ConstraintDifference(*oldconstr);
break;
}
case P2PDistance: {
ConstraintP2PDistance *oldconstr = static_cast<ConstraintP2PDistance *>(*constr);
newconstr = new ConstraintP2PDistance(*oldconstr);
break;
}
case P2PAngle: {
ConstraintP2PAngle *oldconstr = static_cast<ConstraintP2PAngle *>(*constr);
newconstr = new ConstraintP2PAngle(*oldconstr);
break;
}
case P2LDistance: {
ConstraintP2LDistance *oldconstr = static_cast<ConstraintP2LDistance *>(*constr);
newconstr = new ConstraintP2LDistance(*oldconstr);
break;
}
case PointOnLine: {
ConstraintPointOnLine *oldconstr = static_cast<ConstraintPointOnLine *>(*constr);
newconstr = new ConstraintPointOnLine(*oldconstr);
break;
}
case Parallel: {
ConstraintParallel *oldconstr = static_cast<ConstraintParallel *>(*constr);
newconstr = new ConstraintParallel(*oldconstr);
break;
}
case Perpendicular: {
ConstraintPerpendicular *oldconstr = static_cast<ConstraintPerpendicular
*>(*constr); newconstr = new ConstraintPerpendicular(*oldconstr); break;
}
case L2LAngle: {
ConstraintL2LAngle *oldconstr = static_cast<ConstraintL2LAngle *>(*constr);
newconstr = new ConstraintL2LAngle(*oldconstr);
break;
}
case MidpointOnLine: {
ConstraintMidpointOnLine *oldconstr = static_cast<ConstraintMidpointOnLine
*>(*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<Eigen::Infinity>() + 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);

21
src/Mod/Sketcher/App/planegcs/Geo.cpp
View File

@@ -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;
}

7
src/Mod/Sketcher/App/planegcs/Geo.h
View File

@@ -23,10 +23,8 @@
#ifndef PLANEGCS_GEO_H
#define PLANEGCS_GEO_H
#include <cmath>
#include "Util.h"
#include <boost/math/constants/constants.hpp>
namespace GCS
{
@@ -50,6 +48,9 @@ public:
};
using VEC_P = std::vector<Point>;
static constexpr double pi = boost::math::constants::pi<double>();
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

17
src/Mod/Sketcher/App/planegcs/SubSystem.cpp
View File

@@ -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<Constraint *> constrs=p2c[param->second];
for (std::vector<Constraint *>::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<Constraint*> constrs = p2c[pmapfind->second];
for (std::vector<Constraint*>::const_iterator constr = constrs.begin();
constr != constrs.end();

5
src/Mod/Sketcher/App/planegcs/Util.h
View File

@@ -38,11 +38,6 @@ using MAP_pD_D = std::map<double*, double>;
using MAP_pD_I = std::map<double*, int>;
using SET_pD = std::set<double*>;
using SET_I = std::set<int>;
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
} // namespace GCS
#endif // PLANEGCS_UTIL_H