diff --git a/include/ceres/fpclassify.h b/include/ceres/fpclassify.h new file mode 100644 index 000000000..59424bf3e --- /dev/null +++ b/include/ceres/fpclassify.h @@ -0,0 +1,49 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2012 Google Inc. All rights reserved. +// http://code.google.com/p/ceres-solver/ +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// * Neither the name of Google Inc. nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// +// Author: keir@google.com (Keir Mierle) +// +// Portable floating point classification. The names are picked such that they +// do not collide with macros. For example, "isnan" in C99 is a macro and hence +// does not respect namespaces. +// +// TODO(keir): Finish porting! + +#ifndef CERES_PUBLIC_FPCLASSIFY_H_ +#define CERES_PUBLIC_FPCLASSIFY_H_ + +namespace ceres { + +inline bool IsFinite (double x) { return std::isfinite(x); } +inline bool IsInfinite(double x) { return std::isinf(x); } +inline bool IsNaN (double x) { return std::isnan(x); } +inline bool IsNormal (double x) { return std::isnormal(x); } + +} // namespace ceres + +#endif // CERES_PUBLIC_FPCLASSIFY_H_ diff --git a/include/ceres/jet.h b/include/ceres/jet.h index 50b5555c5..f9997bdf2 100644 --- a/include/ceres/jet.h +++ b/include/ceres/jet.h @@ -162,6 +162,7 @@ #include #include "Eigen/Core" +#include "ceres/fpclassify.h" namespace ceres { @@ -401,10 +402,6 @@ inline double cos (double x) { return std::cos(x); } inline double acos (double x) { return std::acos(x); } inline double sin (double x) { return std::sin(x); } inline double asin (double x) { return std::asin(x); } -inline bool isfinite(double x) { return std::isfinite(x); } -inline bool isinf (double x) { return std::isinf(x); } -inline bool isnan (double x) { return std::isnan(x); } -inline bool isnormal(double x) { return std::isnormal(x); } inline double pow (double x, double y) { return std::pow(x, y); } inline double atan2(double y, double x) { return std::atan2(y, x); } @@ -482,22 +479,23 @@ Jet asin(const Jet& f) { } // Jet Classification. It is not clear what the appropriate semantics are for -// these classifications. This picks that isfinite and isnormal are "all" -// operations, i.e. all elements of the jet must be finite for the jet itself to -// be finite (or normal). For isnan and isinf, the answer is less clear. This -// takes a "any" approach for isnan and isinf such that if any part of a jet is -// nan or inf, then the entire jet is nan or inf. This leads to strange -// situations like a jet can be both isinf and isnan, but in practice the "any" -// semantics are the most useful for e.g. checking that derivatives are sane. +// these classifications. This picks that IsFinite and isnormal are "all" +// operations, i.e. all elements of the jet must be finite for the jet itself +// to be finite (or normal). For IsNaN and IsInfinite, the answer is less +// clear. This takes a "any" approach for IsNaN and IsInfinite such that if any +// part of a jet is nan or inf, then the entire jet is nan or inf. This leads +// to strange situations like a jet can be both IsInfinite and IsNaN, but in +// practice the "any" semantics are the most useful for e.g. checking that +// derivatives are sane. // The jet is finite if all parts of the jet are finite. template inline -bool isfinite(const Jet& f) { - if (!isfinite(f.a)) { +bool IsFinite(const Jet& f) { + if (!IsFinite(f.a)) { return false; } for (int i = 0; i < N; ++i) { - if (!isfinite(f.v[i])) { + if (!IsFinite(f.v[i])) { return false; } } @@ -506,12 +504,12 @@ bool isfinite(const Jet& f) { // The jet is infinite if any part of the jet is infinite. template inline -bool isinf(const Jet& f) { - if (isinf(f.a)) { +bool IsInfinite(const Jet& f) { + if (IsInfinite(f.a)) { return true; } for (int i = 0; i < N; i++) { - if (isinf(f.v[i])) { + if (IsFinite(f.v[i])) { return true; } } @@ -520,12 +518,12 @@ bool isinf(const Jet& f) { // The jet is NaN if any part of the jet is NaN. template inline -bool isnan(const Jet& f) { - if (isnan(f.a)) { +bool IsNaN(const Jet& f) { + if (IsNaN(f.a)) { return true; } for (int i = 0; i < N; ++i) { - if (isnan(f.v[i])) { + if (IsNaN(f.v[i])) { return true; } } @@ -534,12 +532,12 @@ bool isnan(const Jet& f) { // The jet is normal if all parts of the jet are normal. template inline -bool isnormal(const Jet& f) { - if (!isnormal(f.a)) { +bool IsNormal(const Jet& f) { + if (!IsNormal(f.a)) { return false; } for (int i = 0; i < N; ++i) { - if (!isnormal(f.v[i])) { + if (!IsNormal(f.v[i])) { return false; } } diff --git a/internal/ceres/array_utils.cc b/internal/ceres/array_utils.cc index aab0c21ed..673baa4f7 100644 --- a/internal/ceres/array_utils.cc +++ b/internal/ceres/array_utils.cc @@ -32,6 +32,7 @@ #include #include +#include "ceres/fpclassify.h" namespace ceres { namespace internal { @@ -46,7 +47,7 @@ const double kImpossibleValue = 1e302; bool IsArrayValid(const int size, const double* x) { if (x != NULL) { for (int i = 0; i < size; ++i) { - if (!isfinite(x[i]) || (x[i] == kImpossibleValue)) { + if (!IsFinite(x[i]) || (x[i] == kImpossibleValue)) { return false; } } diff --git a/internal/ceres/conjugate_gradients_solver.cc b/internal/ceres/conjugate_gradients_solver.cc index 09a02796f..7501bb2db 100644 --- a/internal/ceres/conjugate_gradients_solver.cc +++ b/internal/ceres/conjugate_gradients_solver.cc @@ -42,6 +42,7 @@ #include #include #include +#include "ceres/fpclassify.h" #include "ceres/linear_operator.h" #include "ceres/internal/eigen.h" #include "ceres/types.h" @@ -51,7 +52,7 @@ namespace internal { namespace { bool IsZeroOrInfinity(double x) { - return ((x == 0.0) || (isinf(x))); + return ((x == 0.0) || (IsInfinite(x))); } // Constant used in the MATLAB implementation ~ 2 * eps. @@ -150,14 +151,14 @@ LinearSolver::Summary ConjugateGradientsSolver::Solve( A->RightMultiply(p.data(), q.data()); double pq = p.dot(q); - if ((pq <= 0) || isinf(pq)) { + if ((pq <= 0) || IsInfinite(pq)) { LOG(ERROR) << "Numerical failure. pq = " << pq; summary.termination_type = FAILURE; break; } double alpha = rho / pq; - if (isinf(alpha)) { + if (IsInfinite(alpha)) { LOG(ERROR) << "Numerical failure. alpha " << alpha; summary.termination_type = FAILURE; break; diff --git a/internal/ceres/jet_test.cc b/internal/ceres/jet_test.cc index 2e47b63eb..f12b8b988 100644 --- a/internal/ceres/jet_test.cc +++ b/internal/ceres/jet_test.cc @@ -35,6 +35,7 @@ #include #include "gtest/gtest.h" +#include "ceres/fpclassify.h" #include "ceres/stringprintf.h" #include "ceres/test_util.h" @@ -289,10 +290,10 @@ TEST(JetTraitsTest, ClassificationMixed) { a.v[0] = std::numeric_limits::quiet_NaN(); a.v[1] = std::numeric_limits::infinity(); a.v[2] = -std::numeric_limits::infinity(); - EXPECT_FALSE(isfinite(a)); - EXPECT_FALSE(isnormal(a)); - EXPECT_TRUE(isinf(a)); - EXPECT_TRUE(isnan(a)); + EXPECT_FALSE(IsFinite(a)); + EXPECT_FALSE(IsNormal(a)); + EXPECT_TRUE(IsInfinite(a)); + EXPECT_TRUE(IsNaN(a)); } TEST(JetTraitsTest, ClassificationNaN) { @@ -300,10 +301,10 @@ TEST(JetTraitsTest, ClassificationNaN) { a.v[0] = std::numeric_limits::quiet_NaN(); a.v[1] = 0.0; a.v[2] = 0.0; - EXPECT_FALSE(isfinite(a)); - EXPECT_FALSE(isnormal(a)); - EXPECT_FALSE(isinf(a)); - EXPECT_TRUE(isnan(a)); + EXPECT_FALSE(IsFinite(a)); + EXPECT_FALSE(IsNormal(a)); + EXPECT_FALSE(IsInfinite(a)); + EXPECT_TRUE(IsNaN(a)); } TEST(JetTraitsTest, ClassificationInf) { @@ -311,10 +312,10 @@ TEST(JetTraitsTest, ClassificationInf) { a.v[0] = std::numeric_limits::infinity(); a.v[1] = 0.0; a.v[2] = 0.0; - EXPECT_FALSE(isfinite(a)); - EXPECT_FALSE(isnormal(a)); - EXPECT_TRUE(isinf(a)); - EXPECT_FALSE(isnan(a)); + EXPECT_FALSE(IsFinite(a)); + EXPECT_FALSE(IsNormal(a)); + EXPECT_TRUE(IsInfinite(a)); + EXPECT_FALSE(IsNaN(a)); } TEST(JetTraitsTest, ClassificationFinite) { @@ -322,10 +323,10 @@ TEST(JetTraitsTest, ClassificationFinite) { a.v[0] = 100.0; a.v[1] = 1.0; a.v[2] = 3.14159; - EXPECT_TRUE(isfinite(a)); - EXPECT_TRUE(isnormal(a)); - EXPECT_FALSE(isinf(a)); - EXPECT_FALSE(isnan(a)); + EXPECT_TRUE(IsFinite(a)); + EXPECT_TRUE(IsNormal(a)); + EXPECT_FALSE(IsInfinite(a)); + EXPECT_FALSE(IsNaN(a)); } } // namespace internal diff --git a/internal/ceres/local_parameterization_test.cc b/internal/ceres/local_parameterization_test.cc index 8315e8e2f..aa49250d9 100644 --- a/internal/ceres/local_parameterization_test.cc +++ b/internal/ceres/local_parameterization_test.cc @@ -30,6 +30,7 @@ #include #include "gtest/gtest.h" +#include "ceres/fpclassify.h" #include "ceres/internal/autodiff.h" #include "ceres/internal/eigen.h" #include "ceres/local_parameterization.h" @@ -183,7 +184,7 @@ void QuaternionParameterizationTestHelper(const double* x, double jacobian[12]; param.ComputeJacobian(x, jacobian); for (int i = 0; i < 12; ++i) { - EXPECT_TRUE(isfinite(jacobian[i])); + EXPECT_TRUE(IsFinite(jacobian[i])); EXPECT_NEAR(jacobian[i], jacobian_ref[i], kTolerance) << "Jacobian mismatch: i = " << i << "\n Expected \n" << ConstMatrixRef(jacobian_ref, 4, 3) diff --git a/internal/ceres/rotation_test.cc b/internal/ceres/rotation_test.cc index a8351d50b..55a4fbb99 100644 --- a/internal/ceres/rotation_test.cc +++ b/internal/ceres/rotation_test.cc @@ -591,8 +591,8 @@ J4 MakeJ4(double a, double v0, double v1, double v2, double v3) { bool IsClose(double x, double y) { - EXPECT_FALSE(isnan(x)); - EXPECT_FALSE(isnan(y)); + EXPECT_FALSE(IsNaN(x)); + EXPECT_FALSE(IsNaN(y)); double absdiff = fabs(x - y); if (x == 0 || y == 0) { return absdiff <= kTolerance;