LSSTApplications  10.0+286,10.0+36,10.0+46,10.0-2-g4f67435,10.1+152,10.1+37,11.0,11.0+1,11.0-1-g47edd16,11.0-1-g60db491,11.0-1-g7418c06,11.0-2-g04d2804,11.0-2-g68503cd,11.0-2-g818369d,11.0-2-gb8b8ce7
LSSTDataManagementBasePackage
Public Member Functions | Private Types | Private Member Functions | Private Attributes | List of all members
lsst::meas::algorithms::shapelet::NLSolver Class Referenceabstract

#include <NLSolver.h>

Inheritance diagram for lsst::meas::algorithms::shapelet::NLSolver:
lsst::meas::algorithms::shapelet::CrudeSolver lsst::meas::algorithms::shapelet::EllipseSolver3

Public Member Functions

 NLSolver ()
 
virtual ~NLSolver ()
 
virtual void calculateF (const DVector &x, DVector &f) const =0
 
virtual void calculateJ (const DVector &x, const DVector &f, DMatrix &j) const
 
virtual bool solve (DVector &x, DVector &f) const
 
virtual void getCovariance (DMatrix &cov) const
 
virtual void getInverseCovariance (DMatrix &invcov) const
 
virtual bool testJ (const DVector &, DVector &, std::ostream *os=0, double relerr=0.) const
 
virtual void useHybrid ()
 
virtual void useDogleg ()
 
virtual void setFTol (double fTol)
 
virtual void setGTol (double gTol)
 
virtual void setTol (double fTol, double gTol)
 
virtual void setMinStep (double minStep)
 
virtual void setMaxIter (int maxIter)
 
virtual void setTau (double tau)
 
virtual void setDelta0 (double delta0)
 
virtual double getFTol ()
 
virtual double getGTol ()
 
virtual double getMinStep ()
 
virtual int getMaxIter ()
 
virtual double getTau ()
 
virtual double getDelta0 ()
 
virtual void setOutput (std::ostream &os)
 
virtual void useVerboseOutput ()
 
virtual void useExtraVerboseOutput ()
 
virtual void noUseVerboseOutput ()
 
virtual void useDirectH ()
 
virtual void useSVD ()
 
virtual void useCholesky ()
 
virtual void noUseDirectH ()
 
virtual void noUseSVD ()
 
virtual void noUseCholesky ()
 

Private Types

enum  Method { HYBRID, DOGLEG }
 

Private Member Functions

bool solveDogleg (DVector &x, DVector &f) const
 
bool solveHybrid (DVector &x, DVector &f) const
 

Private Attributes

Method _method
 
double _fTol
 
double _gTol
 
double _minStep
 
int _maxIter
 
double _tau
 
double _delta0
 
std::ostream * _nlOut
 
int _verbose
 
bool _hasDirectH
 
bool _shouldUseCh
 
bool _shouldUseSvd
 
std::auto_ptr< DMatrix > _pJ
 

Detailed Description

Definition at line 356 of file NLSolver.h.

Member Enumeration Documentation

enum lsst::meas::algorithms::shapelet::NLSolver::Method
private
Enumerator
HYBRID 
DOGLEG 

Definition at line 424 of file NLSolver.h.

424 { HYBRID, DOGLEG };

Constructor & Destructor Documentation

lsst::meas::algorithms::shapelet::NLSolver::NLSolver ( )

Definition at line 1194 of file NLSolver.cc.

1194  :
1195  _method(HYBRID),
1196  _fTol(1.e-8), _gTol(1.e-8), _minStep(1.e-8), _maxIter(200),
1197  _tau(1.e-3), _delta0(1.),
1198  _nlOut(0), _verbose(0)
1199  {}
virtual lsst::meas::algorithms::shapelet::NLSolver::~NLSolver ( )
inlinevirtual

Definition at line 361 of file NLSolver.h.

361 {}

Member Function Documentation

virtual void lsst::meas::algorithms::shapelet::NLSolver::calculateF ( const DVector &  x,
DVector &  f 
) const
pure virtual

Implemented in lsst::meas::algorithms::shapelet::CrudeSolver, and lsst::meas::algorithms::shapelet::EllipseSolver3.

void lsst::meas::algorithms::shapelet::NLSolver::calculateJ ( const DVector &  x,
const DVector &  f,
DMatrix &  j 
) const
virtual

Reimplemented in lsst::meas::algorithms::shapelet::CrudeSolver, and lsst::meas::algorithms::shapelet::EllipseSolver3.

Definition at line 1201 of file NLSolver.cc.

1203  {
1204  const double sqrtEps = sqrt(std::numeric_limits<double>::epsilon());
1205  // Do a finite difference calculation for J.
1206  // This function is virtual, so if there is a better way to
1207  // calculate J, then you should override this version.
1208 
1209  DVector x2 = x;
1210  DVector f2(f.size());
1211  DVector f1(f.size());
1212  const int n = x.size();
1213  for(int j=0;j<n;++j) {
1214  const double dx = sqrtEps * (x.norm() + 1.);
1215  x2(j) += dx;
1216  this->calculateF(x2,f2);
1217  x2(j) -= 2.*dx;
1218  this->calculateF(x2,f1);
1219  df.col(j) = (f2-f1)/(2.*dx);
1220  x2(j) = x(j);
1221  }
1222  }
lsst::meas::algorithms::shapelet::NLSolver::calculateF
virtual void calculateF(const DVector &x, DVector &f) const =0
x
double x
Definition: ChebyshevBoundedField.cc:305
void lsst::meas::algorithms::shapelet::NLSolver::getCovariance ( DMatrix &  cov) const
virtual

Reimplemented in lsst::meas::algorithms::shapelet::EllipseSolver3.

Definition at line 1638 of file NLSolver.cc.

1639  {
1640  const double eps = std::numeric_limits<double>::epsilon();
1641  if (!_pJ.get()) {
1642  throw std::runtime_error(
1643  "J not set before calling getCovariance");
1644  }
1645  DMatrix& J = *_pJ;
1646  // This might have changed between solve and getCovariance:
1647  // And we need to set the threshold to sqrt(eps) rather than eps
1648  if (_shouldUseSvd) {
1649  Eigen::JacobiSVD<DMatrix> SV_Solver_J(J, Eigen::ComputeThinU | Eigen::ComputeThinV);
1650  const DVector& svd_s = SV_Solver_J.singularValues();
1651  int kmax = svd_s.size();
1652  while(svd_s(kmax-1) < eps * svd_s(0)) --kmax;
1653  dbg<<"In NLSolver::getCovariance:\n";
1654  dbg<<"Using kmax = "<<kmax<<" size = "<<svd_s.size()<<std::endl;
1655  // (JtJ)^-1 = ( (USVt)t (USVt) )^-1
1656  // = ( V St Ut U S Vt )^-1
1657  // = ( V S^2 Vt )^-1
1658  // = V S^-2 Vt
1659  const DMatrix& svd_v = SV_Solver_J.matrixV();
1660  DVector sm2 = svd_s.array().square().inverse().matrix();
1661  sm2.TMV_subVector(kmax, svd_s.size()).setZero();
1662  cov = svd_v * sm2.asDiagonal() * svd_v.transpose();
1663  } else {
1664  Eigen::ColPivHouseholderQR<DMatrix> QR_Solver_J(J);
1665  // (JtJ)^-1 = ( (QR)t (QR) )^-1
1666  // = ( Rt Qt Q R ) ^-1
1667  // = ( Rt R )^-1
1668  // = R^-1 Rt^-1
1669  cov.setIdentity();
1670  QR_Solver_J.matrixQR().triangularView<Eigen::Upper>().transpose().solveInPlace(cov);
1671  QR_Solver_J.matrixQR().triangularView<Eigen::Upper>().solveInPlace(cov);
1672  }
1673  }
dbg
#define dbg
Definition: NLSolver.cc:37
virtual double lsst::meas::algorithms::shapelet::NLSolver::getDelta0 ( )
inlinevirtual

Definition at line 408 of file NLSolver.h.

408 { return _delta0; }
virtual double lsst::meas::algorithms::shapelet::NLSolver::getFTol ( )
inlinevirtual

Definition at line 403 of file NLSolver.h.

403 { return _fTol; }
virtual double lsst::meas::algorithms::shapelet::NLSolver::getGTol ( )
inlinevirtual

Definition at line 404 of file NLSolver.h.

404 { return _gTol; }
void lsst::meas::algorithms::shapelet::NLSolver::getInverseCovariance ( DMatrix &  invcov) const
virtual

Reimplemented in lsst::meas::algorithms::shapelet::EllipseSolver3.

Definition at line 1675 of file NLSolver.cc.

1676  {
1677  if (!_pJ.get()) {
1678  throw std::runtime_error(
1679  "J not set before calling getInverseCovariance");
1680  }
1681  DMatrix& J = *_pJ;
1682  invCov = J.transpose() * J;
1683  }
virtual int lsst::meas::algorithms::shapelet::NLSolver::getMaxIter ( )
inlinevirtual

Definition at line 406 of file NLSolver.h.

406 { return _maxIter; }
virtual double lsst::meas::algorithms::shapelet::NLSolver::getMinStep ( )
inlinevirtual

Definition at line 405 of file NLSolver.h.

405 { return _minStep; }
virtual double lsst::meas::algorithms::shapelet::NLSolver::getTau ( )
inlinevirtual

Definition at line 407 of file NLSolver.h.

407 { return _tau; }
virtual void lsst::meas::algorithms::shapelet::NLSolver::noUseCholesky ( )
inlinevirtual

Definition at line 420 of file NLSolver.h.

420 {}
virtual void lsst::meas::algorithms::shapelet::NLSolver::noUseDirectH ( )
inlinevirtual

Definition at line 418 of file NLSolver.h.

418 {}
virtual void lsst::meas::algorithms::shapelet::NLSolver::noUseSVD ( )
inlinevirtual

Definition at line 419 of file NLSolver.h.

419 { _shouldUseSvd = false; }
virtual void lsst::meas::algorithms::shapelet::NLSolver::noUseVerboseOutput ( )
inlinevirtual

Definition at line 413 of file NLSolver.h.

413 { _verbose = 0; }
virtual void lsst::meas::algorithms::shapelet::NLSolver::setDelta0 ( double  delta0)
inlinevirtual

Definition at line 401 of file NLSolver.h.

401 { _delta0 = delta0; }
virtual void lsst::meas::algorithms::shapelet::NLSolver::setFTol ( double  fTol)
inlinevirtual

Definition at line 393 of file NLSolver.h.

393 { _fTol = fTol; }
virtual void lsst::meas::algorithms::shapelet::NLSolver::setGTol ( double  gTol)
inlinevirtual

Definition at line 394 of file NLSolver.h.

394 { _gTol = gTol; }
virtual void lsst::meas::algorithms::shapelet::NLSolver::setMaxIter ( int  maxIter)
inlinevirtual

Definition at line 399 of file NLSolver.h.

399 { _maxIter = maxIter; }
virtual void lsst::meas::algorithms::shapelet::NLSolver::setMinStep ( double  minStep)
inlinevirtual

Definition at line 398 of file NLSolver.h.

398 { _minStep = minStep; }
virtual void lsst::meas::algorithms::shapelet::NLSolver::setOutput ( std::ostream &  os)
inlinevirtual

Definition at line 410 of file NLSolver.h.

410 { _nlOut = &os; }
virtual void lsst::meas::algorithms::shapelet::NLSolver::setTau ( double  tau)
inlinevirtual

Definition at line 400 of file NLSolver.h.

400 { _tau = tau; }
virtual void lsst::meas::algorithms::shapelet::NLSolver::setTol ( double  fTol,
double  gTol 
)
inlinevirtual

Definition at line 395 of file NLSolver.h.

396  { _fTol = fTol; _gTol = gTol; }
bool lsst::meas::algorithms::shapelet::NLSolver::solve ( DVector &  x,
DVector &  f 
) const
virtual

Reimplemented in lsst::meas::algorithms::shapelet::EllipseSolver3.

Definition at line 1620 of file NLSolver.cc.

1621  {
1622 #ifndef NOTHROW
1623  try {
1624 #endif
1625  switch (_method) {
1626  case HYBRID : return solveHybrid(x,f);
1627  case DOGLEG : return solveDogleg(x,f);
1628  default : dbg<<"Unknown method\n"; return false;
1629  }
1630 #ifndef NOTHROW
1631  } catch (...) {
1632  dbg<<"Error encountered in NLSolver::Solve\n";
1633  }
1634 #endif
1635  return false;
1636  }
lsst::meas::algorithms::shapelet::NLSolver::solveDogleg
bool solveDogleg(DVector &x, DVector &f) const
Definition: NLSolver.cc:1307
lsst::meas::algorithms::shapelet::NLSolver::solveHybrid
bool solveHybrid(DVector &x, DVector &f) const
Definition: NLSolver.cc:1426
x
double x
Definition: ChebyshevBoundedField.cc:305
dbg
#define dbg
Definition: NLSolver.cc:37
bool lsst::meas::algorithms::shapelet::NLSolver::solveDogleg ( DVector &  x,
DVector &  f 
) const
private

Definition at line 1307 of file NLSolver.cc.

1309  {
1310  dbg<<"Start Solve_Dogleg\n";
1311  _pJ.reset(new DMatrix(f.size(),x.size()));
1312  DMatrix& J = *_pJ;
1313  DVector h(x.size());
1314  DVector temp(x.size());
1315  DVector xNew(x.size());
1316  DVector fNew(f.size());
1317 
1318  xdbg<<"x = "<<x.transpose()<<std::endl;
1319  this->calculateF(x,f);
1320  xdbg<<"f = "<<f.transpose()<<std::endl;
1321  CHECKF(f.TMV_normInf());
1322  double Q = 0.5*f.TMV_normSq();
1323  xdbg<<"Q = "<<Q<<std::endl;
1324  this->calculateJ(x,f,J);
1325  xdbg<<"J = "<<J<<std::endl;
1326 
1327  DVector g = J.transpose() * f;
1328  xdbg<<"g = "<<g.transpose()<<std::endl;
1329  CHECKG(g.TMV_normInf());
1330 
1331  double delta = _delta0;
1332  int maxnsing = std::min(f.size(),x.size());
1333  int nsing = maxnsing;
1334 
1335  dbg<<"iter |f|inf Q |g|inf delta\n";
1336  for(int k=0;k<_maxIter;++k) {
1337  dbg<<k<<" "<<f.TMV_normInf()<<" "<<Q<<" "<<g.TMV_normInf()<<" "<<delta<<std::endl;
1338  //h = -f/J;
1339  h = J.lu().solve(-f);
1340  xdbg<<"h = "<<h.transpose()<<std::endl;
1341 
1342  double normsqg = g.TMV_normSq();
1343  double normH = h.norm();
1344  double rhoDenom;
1345 
1346  if (normH <= delta) {
1347  xdbg<<"|h| < delta\n";
1348  rhoDenom = Q;
1349  xdbg<<"rhodenom = "<<rhoDenom<<std::endl;
1350  } else {
1351  double alpha = normsqg / (J*g).TMV_normSq();
1352  double normG = sqrt(normsqg);
1353  if (normG >= delta / alpha) {
1354  xdbg<<"|g| > delta/alpha\n";
1355  h = -(delta / normG) * g;
1356  xdbg<<"h = "<<h.transpose()<<std::endl;
1357  rhoDenom = delta*(2.*alpha*normG-delta)/(2.*alpha);
1358  xdbg<<"rhodenom = "<<rhoDenom<<std::endl;
1359  } else {
1360  xdbg<<"dogleg\n";
1361  temp = h + alpha*g;
1362  double a = temp.TMV_normSq();
1363  double b = -alpha * (g.transpose() * temp)(0,0);
1364  double c = alpha*alpha*g.TMV_normSq()-delta*delta;
1365  // beta is the solution of 0 = a beta^2 + 2b beta + c
1366  xdbg<<"a,b,c = "<<a<<" "<<b<<" "<<c<<std::endl;
1367  double beta = (b <= 0) ?
1368  (-b + sqrt(b*b - a*c)) / a :
1369  -c / (b + sqrt(b*b - a*c));
1370  xdbg<<"alpha = "<<alpha<<std::endl;
1371  xdbg<<"beta = "<<beta<<std::endl;
1372  h = -alpha*g + beta*temp;
1373  xdbg<<"h = "<<h.transpose()<<std::endl;
1374  xdbg<<"h.norm() = "<<h.norm()<<" delta = "<<delta<<std::endl;
1375  rhoDenom =
1376  0.5*alpha*std::pow((1.-beta)*normG,2) +
1377  beta*(2.-beta)*Q;
1378  xdbg<<"rhodenom = "<<rhoDenom<<std::endl;
1379  }
1380  normH = h.norm();
1381  }
1382 
1383  CHECKSTEP(normH);
1384 
1385  xNew = x + h;
1386  xdbg<<"xnew = "<<xNew.transpose()<<std::endl;
1387  this->calculateF(xNew,fNew);
1388  xdbg<<"fnew = "<<fNew.transpose()<<std::endl;
1389  double QNew = 0.5*fNew.TMV_normSq();
1390  xdbg<<"Qnew = "<<QNew<<std::endl;
1391 
1392  bool deltaok = false;
1393  if (QNew < Q) {
1394  double rho = (Q-QNew) / rhoDenom;
1395  xdbg<<"rho = "<<Q-QNew<<" / "<<rhoDenom<<" = "<<rho<<std::endl;
1396  x = xNew; f = fNew; Q = QNew;
1397  CHECKF(f.TMV_normInf());
1398  this->calculateJ(x,f,J);
1399  g = J.transpose() * f;
1400  xdbg<<"g = "<<g.transpose()<<std::endl;
1401  CHECKG(g.TMV_normInf());
1402  if (rho > 0.75) {
1403  delta = std::max(delta,3.*normH);
1404  deltaok = true;
1405  }
1406  }
1407  if (deltaok) {
1408  nsing = maxnsing;
1409  } else {
1410  delta /= 2.;
1411  double min_delta = _minStep * (x.norm()+_minStep);
1412  if (delta < min_delta) {
1413  dbg<<"delta became too small ("<<
1414  delta<<" < "<<min_delta<<")\n";
1415  SHOWFAILFG;
1416  return false;
1417  }
1418  }
1419  }
1420  assert(nsing >= 0); // make compiler happy about unused variable
1421  dbg<<"Maximum iterations exceeded in Dogleg method\n";
1422  SHOWFAILFG;
1423  return false;
1424  }
CHECKG
#define CHECKG(normInfG)
Definition: NLSolver.cc:1279
lsst::meas::algorithms::shapelet::NLSolver::calculateF
virtual void calculateF(const DVector &x, DVector &f) const =0
xdbg
#define xdbg
Definition: NLSolver.cc:38
alpha
Matrix alpha
Definition: GaussianCentroid.cc:76
CHECKF
#define CHECKF(normInfF)
Definition: NLSolver.cc:1269
TMV_normSq
#define TMV_normSq()
Definition: MyMatrix.h:186
lsst::meas::algorithms::shapelet::NLSolver::calculateJ
virtual void calculateJ(const DVector &x, const DVector &f, DMatrix &j) const
Definition: NLSolver.cc:1201
CHECKSTEP
#define CHECKSTEP(normH)
Definition: NLSolver.cc:1295
x
double x
Definition: ChebyshevBoundedField.cc:305
dbg
#define dbg
Definition: NLSolver.cc:37
SHOWFAILFG
#define SHOWFAILFG
Definition: NLSolver.cc:1289
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:49
beta
Vector beta
Definition: GaussianCentroid.cc:77
bool lsst::meas::algorithms::shapelet::NLSolver::solveHybrid ( DVector &  x,
DVector &  f 
) const
private

Definition at line 1426 of file NLSolver.cc.

1429  {
1430  const double sqrtEps = sqrt(std::numeric_limits<double>::epsilon());
1431 
1432  dbg<<"Start Solve_Hybrid\n";
1433  _pJ.reset(new DMatrix(f.size(),x.size()));
1434  DMatrix& J = *_pJ;
1435  DVector h(x.size());
1436  DVector xNew(x.size());
1437  DVector fNew(f.size());
1438  DVector gNew(x.size());
1439  DMatrix JNew(f.size(),x.size());
1440  DVector y(x.size());
1441  DVector v(x.size());
1442 
1443  xdbg<<"x = "<<x.transpose()<<std::endl;
1444  this->calculateF(x,f);
1445  xdbg<<"f = "<<f.transpose()<<std::endl;
1446  double normInfF = f.TMV_normInf();
1447  CHECKF(normInfF);
1448  double Q = 0.5*f.TMV_normSq();
1449  xdbg<<"Q = "<<Q<<std::endl;
1450  this->calculateJ(x,f,J);
1451  xdbg<<"J = "<<J<<std::endl;
1452 
1453  DMatrix A = J.transpose()*J;
1454  xdbg<<"A = "<<A<<std::endl;
1455  DMatrix H(x.size(),x.size());
1456  xdbg<<"setToIdent\n";
1457  H.TMV_setToIdentity();
1458  xdbg<<"After calculate H = "<<H<<std::endl;
1459 
1460  DVector g = J.transpose() * f;
1461  xdbg<<"g = "<<g.transpose()<<std::endl;
1462  double normInfG = g.TMV_normInf();
1463  CHECKG(normInfG);
1464 
1465  double mu = _tau * A.TMV_diag().TMV_normInf();
1466  A.EIGEN_diag() += mu;
1467  double nu = 2.;
1468  double delta = _delta0;
1469  bool shouldUseQuasiNewton = false;
1470  int count = 0;
1471 
1472  dbg<<"iter |f|inf Q |g|inf mu delta LM/QN\n";
1473  for(int k=0;k<_maxIter;++k) {
1474  dbg<<k<<" "<<normInfF<<" "<<Q<<" "<<normInfG<<" "<< mu<<" "<<
1475  delta<<" "<<(shouldUseQuasiNewton?"QN":"LM")<<std::endl;
1476  xdbg<<"k = "<<k<<std::endl;
1477  xdbg<<"mu = "<<mu<<std::endl;
1478  xdbg<<"delta = "<<delta<<std::endl;
1479  xdbg<<"A = "<<A<<std::endl;
1480  xdbg<<"H = "<<H<<std::endl;
1481  xdbg<<"method = "<<(shouldUseQuasiNewton ? "quasinewton\n" : "LM\n");
1482  bool isBetter = false;
1483  bool shouldSwitchMethod = false;
1484 
1485  if (shouldUseQuasiNewton) {
1486  //h = -g/H;
1487  h = H.ldlt().solve(-g);
1488  } else {
1489  //h = -g/A;
1490  h = A.ldlt().solve(-g);
1491  }
1492 
1493  xdbg<<"h = "<<h.transpose()<<std::endl;
1494  double normH = h.norm();
1495  CHECKSTEP(normH);
1496  if (shouldUseQuasiNewton && normH > delta) h *= delta/normH;
1497 
1498  xNew = x + h;
1499  xdbg<<"xnew = "<<xNew.transpose()<<std::endl;
1500  this->calculateF(xNew,fNew);
1501  xdbg<<"fnew = "<<fNew.transpose()<<std::endl;
1502  double QNew = 0.5*fNew.TMV_normSq();
1503  xdbg<<"Qnew = "<<QNew<<std::endl;
1504 
1505  this->calculateJ(xNew,fNew,JNew);
1506  xdbg<<"Jnew = "<<JNew<<std::endl;
1507  if (shouldUseQuasiNewton || QNew < Q) {
1508  gNew = JNew.transpose() * fNew;
1509  xdbg<<"gnew = "<<gNew.transpose()<<std::endl;
1510  }
1511  double normInfGNew = gNew.TMV_normInf();
1512 
1513  if (shouldUseQuasiNewton) {
1514  xdbg<<"quasinewton\n";
1515  isBetter =
1516  (QNew < Q) ||
1517  (QNew <= (1.+sqrtEps)*Q && normInfGNew < normInfG);
1518  xdbg<<"better = "<<isBetter<<std::endl;
1519  shouldSwitchMethod = (normInfGNew >= normInfG);
1520  xdbg<<"switchmethod = "<<shouldSwitchMethod<<std::endl;
1521  if (QNew < Q) {
1522  double rho = (Q-QNew) / (-(h.transpose()*g)(0,0)-0.5*(J*h).TMV_normSq());
1523  if (rho > 0.75) {
1524  delta = std::max(delta,3.*normH);
1525  } else if (rho < 0.25) {
1526  delta /= 2.;
1527  double min_delta = _minStep * (x.norm()+_minStep);
1528  if (delta < min_delta) {
1529  dbg<<"delta became too small ("<<
1530  delta<<" < "<<min_delta<<")\n";
1531  SHOWFAILFG;
1532  return false;
1533  }
1534  }
1535  } else {
1536  delta /= 2.;
1537  double min_delta = _minStep * (x.norm()+_minStep);
1538  if (delta < min_delta) {
1539  dbg<<"delta became too small ("<<
1540  delta<<" < "<<min_delta<<")\n";
1541  SHOWFAILFG;
1542  return false;
1543  }
1544  }
1545  } else {
1546  xdbg<<"LM\n";
1547  if (QNew <= Q) {
1548  isBetter = true;
1549  // we don't need the g vector anymore, so use this space
1550  // to calculate g-mu*h
1551  //double rho = (Q-QNew) / (0.5*h*(mu*h-g));
1552  g -= mu*h;
1553  double rho = (Q-QNew) / (-0.5*(h.transpose()*g)(0,0));
1554  mu *= std::max(1./3.,1.-std::pow(2.*rho-1.,3)); nu = 2.;
1555  xdbg<<"check1: "<<normInfGNew<<" <? "<<0.02*QNew<<std::endl;
1556  xdbg<<"check2: "<<Q-QNew<<" <? "<<0.02*QNew<<std::endl;
1557  if (std::min(normInfGNew,Q-QNew) < 0.02 * QNew) {
1558  ++count;
1559  if (count == 3) shouldSwitchMethod = true;
1560  } else {
1561  count = 0;
1562  }
1563  if (count != 3) {
1564  A = JNew.transpose() * JNew;
1565  A.EIGEN_diag() += mu;
1566  }
1567  } else {
1568  A.EIGEN_diag() += mu*(nu-1.); mu *= nu; nu *= 2.;
1569  count = 0;
1570  // MJ: try this?
1571  shouldSwitchMethod = (nu >= 32.);
1572  }
1573  xdbg<<"better = "<<isBetter<<std::endl;
1574  xdbg<<"switchmethod = "<<shouldSwitchMethod<<std::endl;
1575  xdbg<<"count = "<<count<<std::endl;
1576  }
1577 
1578  y = JNew.transpose()*(JNew*h) + (JNew-J).transpose()*fNew;
1579  double hy = (h.transpose()*y)(0,0);
1580  xdbg<<"hy = "<<hy<<std::endl;
1581  if (hy > 0.) {
1582  v = H*h;
1583  xdbg<<"v = "<<v.transpose()<<std::endl;
1584  xdbg<<"y = "<<y.transpose()<<std::endl;
1585  double hv = (h.transpose()*v)(0,0);
1586  xdbg<<"hv = "<<hv<<std::endl;
1587  H -= (1./hv) * (v * v.transpose());
1588  xdbg<<"H -> "<<H<<std::endl;
1589  H += (1./hy) * (y * y.transpose());
1590  xdbg<<"H -> "<<H<<std::endl;
1591  }
1592 
1593  if (isBetter) {
1594  xdbg<<"better"<<std::endl;
1595  x = xNew; f = fNew; Q = QNew; J = JNew; g = gNew;
1596  normInfF = f.TMV_normInf(); normInfG = normInfGNew;
1597  CHECKF(normInfF);
1598  CHECKG(normInfG);
1599  }
1600  if (shouldSwitchMethod) {
1601  if (shouldUseQuasiNewton) {
1602  xdbg<<"switch to LM\n";
1603  A = J.transpose() * J;
1604  //mu = _tau * A.diag().normInf();
1605  A.EIGEN_diag() += mu;
1606  shouldUseQuasiNewton = false;
1607  count = 0;
1608  } else {
1609  xdbg<<"switch to quasinewton\n";
1610  delta = std::max(1.5*_minStep*(x.norm()+_minStep),0.2*normH);
1611  shouldUseQuasiNewton = true;
1612  }
1613  }
1614  }
1615  dbg<<"Maximum iterations exceeded in Hybrid method\n";
1616  SHOWFAILFG;
1617  return false;
1618  }
y
int y
Definition: GaussianCentroid.cc:36
CHECKG
#define CHECKG(normInfG)
Definition: NLSolver.cc:1279
lsst::meas::algorithms::shapelet::NLSolver::calculateF
virtual void calculateF(const DVector &x, DVector &f) const =0
xdbg
#define xdbg
Definition: NLSolver.cc:38
CHECKF
#define CHECKF(normInfF)
Definition: NLSolver.cc:1269
TMV_normSq
#define TMV_normSq()
Definition: MyMatrix.h:186
lsst::meas::algorithms::shapelet::NLSolver::calculateJ
virtual void calculateJ(const DVector &x, const DVector &f, DMatrix &j) const
Definition: NLSolver.cc:1201
CHECKSTEP
#define CHECKSTEP(normH)
Definition: NLSolver.cc:1295
x
double x
Definition: ChebyshevBoundedField.cc:305
dbg
#define dbg
Definition: NLSolver.cc:37
SHOWFAILFG
#define SHOWFAILFG
Definition: NLSolver.cc:1289
bool lsst::meas::algorithms::shapelet::NLSolver::testJ ( const DVector &  x,
DVector &  f,
std::ostream *  os = 0,
double  relerr = 0. 
) const
virtual

Reimplemented in lsst::meas::algorithms::shapelet::EllipseSolver3.

Definition at line 1224 of file NLSolver.cc.

1227  {
1228  const double sqrtEps = sqrt(std::numeric_limits<double>::epsilon());
1229 
1230  this->calculateF(x,f);
1231  _pJ.reset(new DMatrix(f.size(),x.size()));
1232  DMatrix& J = *_pJ;
1233  this->calculateJ(x,f,J);
1234  DMatrix Jn(f.size(),x.size());
1235  NLSolver::calculateJ(x,f,Jn);
1236  double err = (J-Jn).TMV_maxAbsElement() / Jn.norm();
1237  if (!relErr) relErr = 10.*sqrtEps;
1238  if (os) {
1239  *os << "TestJ:\n";
1240  if (_verbose >= 1) {
1241  *os << "x = "<<x<<std::endl;
1242  *os << "f = "<<f<<std::endl;
1243  *os << "Direct J = "<<J<<std::endl;
1244  *os << "Numeric J = "<<Jn<<std::endl;
1245  }
1246  *os << "MaxAbsElement(J-J_num) / J.norm() = "<<err<<std::endl;
1247  *os << "cf. relerr = "<<relErr<<std::endl;
1248  if (err >= relErr) {
1249  DMatrix diff = J-Jn;
1250  *os << "J-J_num = "<<diff;
1251  double maxel = diff.TMV_maxAbsElement();
1252  *os << "Max element = "<<maxel<<std::endl;
1253  const int m = diff.TMV_colsize();
1254  const int n = diff.TMV_rowsize();
1255  for(int i=0;i<m;++i) {
1256  for(int j=0;j<n;++j) {
1257  if (std::abs(diff(i,j)) > 0.9*maxel) {
1258  *os<<"J("<<i<<','<<j<<") = "<<J(i,j)<<" ";
1259  *os<<"J_num("<<i<<','<<j<<") = "<<Jn(i,j)<<" ";
1260  *os<<"diff = "<<J(i,j)-Jn(i,j)<<std::endl;
1261  }
1262  }
1263  }
1264  }
1265  }
1266  return err < relErr;
1267  }
lsst::meas::algorithms::shapelet::NLSolver::calculateF
virtual void calculateF(const DVector &x, DVector &f) const =0
lsst::meas::algorithms::shapelet::NLSolver::calculateJ
virtual void calculateJ(const DVector &x, const DVector &f, DMatrix &j) const
Definition: NLSolver.cc:1201
x
double x
Definition: ChebyshevBoundedField.cc:305
TMV_maxAbsElement
#define TMV_maxAbsElement()
Definition: MyMatrix.h:188
lsstimport.m
tuple m
Definition: lsstimport.py:48
virtual void lsst::meas::algorithms::shapelet::NLSolver::useCholesky ( )
inlinevirtual

Definition at line 417 of file NLSolver.h.

417 {}
virtual void lsst::meas::algorithms::shapelet::NLSolver::useDirectH ( )
inlinevirtual

Definition at line 415 of file NLSolver.h.

415 {}
virtual void lsst::meas::algorithms::shapelet::NLSolver::useDogleg ( )
inlinevirtual

Definition at line 391 of file NLSolver.h.

391 { _method = DOGLEG; }
virtual void lsst::meas::algorithms::shapelet::NLSolver::useExtraVerboseOutput ( )
inlinevirtual

Definition at line 412 of file NLSolver.h.

412 { _verbose = 2; }
virtual void lsst::meas::algorithms::shapelet::NLSolver::useHybrid ( )
inlinevirtual

Definition at line 390 of file NLSolver.h.

390 { _method = HYBRID; }
virtual void lsst::meas::algorithms::shapelet::NLSolver::useSVD ( )
inlinevirtual

Definition at line 416 of file NLSolver.h.

416 { _shouldUseSvd = true; }
virtual void lsst::meas::algorithms::shapelet::NLSolver::useVerboseOutput ( )
inlinevirtual

Definition at line 411 of file NLSolver.h.

411 { _verbose = 1; }

Member Data Documentation

double lsst::meas::algorithms::shapelet::NLSolver::_delta0
private

Definition at line 432 of file NLSolver.h.

double lsst::meas::algorithms::shapelet::NLSolver::_fTol
private

Definition at line 427 of file NLSolver.h.

double lsst::meas::algorithms::shapelet::NLSolver::_gTol
private

Definition at line 428 of file NLSolver.h.

bool lsst::meas::algorithms::shapelet::NLSolver::_hasDirectH
private

Definition at line 435 of file NLSolver.h.

int lsst::meas::algorithms::shapelet::NLSolver::_maxIter
private

Definition at line 430 of file NLSolver.h.

Method lsst::meas::algorithms::shapelet::NLSolver::_method
private

Definition at line 426 of file NLSolver.h.

double lsst::meas::algorithms::shapelet::NLSolver::_minStep
private

Definition at line 429 of file NLSolver.h.

std::ostream* lsst::meas::algorithms::shapelet::NLSolver::_nlOut
private

Definition at line 433 of file NLSolver.h.

std::auto_ptr<DMatrix> lsst::meas::algorithms::shapelet::NLSolver::_pJ
mutableprivate

Definition at line 442 of file NLSolver.h.

bool lsst::meas::algorithms::shapelet::NLSolver::_shouldUseCh
private

Definition at line 436 of file NLSolver.h.

bool lsst::meas::algorithms::shapelet::NLSolver::_shouldUseSvd
private

Definition at line 437 of file NLSolver.h.

double lsst::meas::algorithms::shapelet::NLSolver::_tau
private

Definition at line 431 of file NLSolver.h.

int lsst::meas::algorithms::shapelet::NLSolver::_verbose
private

Definition at line 434 of file NLSolver.h.

The documentation for this class was generated from the following files:
Generated on Thu Sep 24 2015 02:29:38 for LSSTApplications by   doxygen 1.8.5