Home | History | Annotate | Download | only in ceres
      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 (at) 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 
     38 namespace ceres {
     39 
     40 void TrivialLoss::Evaluate(double s, double rho[3]) const {
     41   rho[0] = s;
     42   rho[1] = 1;
     43   rho[2] = 0;
     44 }
     45 
     46 void 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 * a_ * r - b_;
     52     rho[1] = a_ / r;
     53     rho[2] = - rho[1] / (2 * s);
     54   } else {
     55     // Inlier region.
     56     rho[0] = s;
     57     rho[1] = 1;
     58     rho[2] = 0;
     59   }
     60 }
     61 
     62 void SoftLOneLoss::Evaluate(double s, double rho[3]) const {
     63   const double sum = 1 + s * c_;
     64   const double tmp = sqrt(sum);
     65   // 'sum' and 'tmp' are always positive, assuming that 's' is.
     66   rho[0] = 2 * b_ * (tmp - 1);
     67   rho[1] = 1 / tmp;
     68   rho[2] = - (c_ * rho[1]) / (2 * sum);
     69 }
     70 
     71 void CauchyLoss::Evaluate(double s, double rho[3]) const {
     72   const double sum = 1 + s * c_;
     73   const double inv = 1 / sum;
     74   // 'sum' and 'inv' are always positive, assuming that 's' is.
     75   rho[0] = b_ * log(sum);
     76   rho[1] = inv;
     77   rho[2] = - c_ * (inv * inv);
     78 }
     79 
     80 void 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] = inv;
     86   rho[2] = -2 * s * b_ * (inv * inv);
     87 }
     88 
     89 TolerantLoss::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 
     97 void 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] = e_x / (1.0 + e_x);
    112     rho[2] = 0.5 / (b_ * (1.0 + cosh(x)));
    113   }
    114 }
    115 
    116 ComposedLoss::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 
    124 ComposedLoss::~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 
    133 void 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 
    144 void 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