1// Ceres Solver - A fast non-linear least squares minimizer
2// Copyright 2010, 2011, 2012 Google Inc. All rights reserved.
3// http://code.google.com/p/ceres-solver/
4//
5// Redistribution and use in source and binary forms, with or without
6// modification, are permitted provided that the following conditions are met:
7//
8// * Redistributions of source code must retain the above copyright notice,
9//   this list of conditions and the following disclaimer.
10// * Redistributions in binary form must reproduce the above copyright notice,
11//   this list of conditions and the following disclaimer in the documentation
12//   and/or other materials provided with the distribution.
13// * Neither the name of Google Inc. nor the names of its contributors may be
14//   used to endorse or promote products derived from this software without
15//   specific prior written permission.
16//
17// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
21// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27// POSSIBILITY OF SUCH DAMAGE.
28//
29// Author: sameeragarwal@google.com (Sameer Agarwal)
30//
31// Purpose: See .h file.
32
33#include "ceres/loss_function.h"
34
35#include <cmath>
36#include <cstddef>
37
38namespace ceres {
39
40void TrivialLoss::Evaluate(double s, double rho[3]) const {
41  rho[0] = s;
42  rho[1] = 1.0;
43  rho[2] = 0.0;
44}
45
46void HuberLoss::Evaluate(double s, double rho[3]) const {
47  if (s > b_) {
48    // Outlier region.
49    // 'r' is always positive.
50    const double r = sqrt(s);
51    rho[0] = 2.0 * a_ * r - b_;
52    rho[1] = std::max(std::numeric_limits<double>::min(), a_ / r);
53    rho[2] = - rho[1] / (2.0 * s);
54  } else {
55    // Inlier region.
56    rho[0] = s;
57    rho[1] = 1.0;
58    rho[2] = 0.0;
59  }
60}
61
62void SoftLOneLoss::Evaluate(double s, double rho[3]) const {
63  const double sum = 1.0 + s * c_;
64  const double tmp = sqrt(sum);
65  // 'sum' and 'tmp' are always positive, assuming that 's' is.
66  rho[0] = 2.0 * b_ * (tmp - 1.0);
67  rho[1] = std::max(std::numeric_limits<double>::min(), 1.0 / tmp);
68  rho[2] = - (c_ * rho[1]) / (2.0 * sum);
69}
70
71void CauchyLoss::Evaluate(double s, double rho[3]) const {
72  const double sum = 1.0 + s * c_;
73  const double inv = 1.0 / sum;
74  // 'sum' and 'inv' are always positive, assuming that 's' is.
75  rho[0] = b_ * log(sum);
76  rho[1] = std::max(std::numeric_limits<double>::min(), inv);
77  rho[2] = - c_ * (inv * inv);
78}
79
80void ArctanLoss::Evaluate(double s, double rho[3]) const {
81  const double sum = 1 + s * s * b_;
82  const double inv = 1 / sum;
83  // 'sum' and 'inv' are always positive.
84  rho[0] = a_ * atan2(s, a_);
85  rho[1] = std::max(std::numeric_limits<double>::min(), inv);
86  rho[2] = -2.0 * s * b_ * (inv * inv);
87}
88
89TolerantLoss::TolerantLoss(double a, double b)
90    : a_(a),
91      b_(b),
92      c_(b * log(1.0 + exp(-a / b))) {
93  CHECK_GE(a, 0.0);
94  CHECK_GT(b, 0.0);
95}
96
97void TolerantLoss::Evaluate(double s, double rho[3]) const {
98  const double x = (s - a_) / b_;
99  // The basic equation is rho[0] = b ln(1 + e^x).  However, if e^x is too
100  // large, it will overflow.  Since numerically 1 + e^x == e^x when the
101  // x is greater than about ln(2^53) for doubles, beyond this threshold
102  // we substitute x for ln(1 + e^x) as a numerically equivalent approximation.
103  static const double kLog2Pow53 = 36.7;  // ln(MathLimits<double>::kEpsilon).
104  if (x > kLog2Pow53) {
105    rho[0] = s - a_ - c_;
106    rho[1] = 1.0;
107    rho[2] = 0.0;
108  } else {
109    const double e_x = exp(x);
110    rho[0] = b_ * log(1.0 + e_x) - c_;
111    rho[1] = std::max(std::numeric_limits<double>::min(), e_x / (1.0 + e_x));
112    rho[2] = 0.5 / (b_ * (1.0 + cosh(x)));
113  }
114}
115
116ComposedLoss::ComposedLoss(const LossFunction* f, Ownership ownership_f,
117                           const LossFunction* g, Ownership ownership_g)
118    : f_(CHECK_NOTNULL(f)),
119      g_(CHECK_NOTNULL(g)),
120      ownership_f_(ownership_f),
121      ownership_g_(ownership_g) {
122}
123
124ComposedLoss::~ComposedLoss() {
125  if (ownership_f_ == DO_NOT_TAKE_OWNERSHIP) {
126    f_.release();
127  }
128  if (ownership_g_ == DO_NOT_TAKE_OWNERSHIP) {
129    g_.release();
130  }
131}
132
133void ComposedLoss::Evaluate(double s, double rho[3]) const {
134  double rho_f[3], rho_g[3];
135  g_->Evaluate(s, rho_g);
136  f_->Evaluate(rho_g[0], rho_f);
137  rho[0] = rho_f[0];
138  // f'(g(s)) * g'(s).
139  rho[1] = rho_f[1] * rho_g[1];
140  // f''(g(s)) * g'(s) * g'(s) + f'(g(s)) * g''(s).
141  rho[2] = rho_f[2] * rho_g[1] * rho_g[1] + rho_f[1] * rho_g[2];
142}
143
144void ScaledLoss::Evaluate(double s, double rho[3]) const {
145  if (rho_.get() == NULL) {
146    rho[0] = a_ * s;
147    rho[1] = a_;
148    rho[2] = 0.0;
149  } else {
150    rho_->Evaluate(s, rho);
151    rho[0] *= a_;
152    rho[1] *= a_;
153    rho[2] *= a_;
154  }
155}
156
157}  // namespace ceres
158