cpp_llsq_obj.hpp

View page source

C++ llsq_obj: Source Code

// ----------------------------------------------------------------------------
# include <cassert>
# include <cmpad/fun_obj.hpp>
//
# include <cmpad/vector.hpp>

namespace cmpad { // BEGIN cmpad namespace

// BEGIN CLASS_DECLARE
template <class Vector> class llsq_obj : public fun_obj<Vector>
// END CLASS_DECLARE
{
private:
   //
   // option_
   option_t option_;
   //
   // t_, q_, y_
   Vector t_, q_, y_;
   //
public:
   // scalar_type
   typedef typename Vector::value_type scalar_type;
   //
   // vector_type
   typedef Vector vector_type;
   //
   // option
   const option_t& option(void) const override
   {  return option_; }
   //
   // domain
   size_t domain(void) const override
   {  return option_.n_arg; }
   //
   // range
   size_t range(void) const override
   {  return 1; }
   //
   // setup
   void setup(const option_t& option) override
   {  //
      // n_arg, n_other
      assert( option.n_arg > 0 );
      assert( option.n_other > 0 );
      //
      // option_
      option_ = option;
      //
      // n_other
      size_t n_other = option.n_other;
      //
      // y_
      y_.resize(1);
      //
      // t_
      t_.resize(n_other);
      if( n_other == 1 )
         t_[0] = 0.0;
      else
      {  for(size_t j = 0; j < n_other; ++j)
            t_[j] = -1.0 + 2.0 * double(j) / double(n_other-1);
      }
      //
      // q_
      q_.resize(n_other);
      if( n_other == 1 )
         q_[0] = 0.0;
      else
      {  for(size_t j = 0; j < n_other; ++j)
            if( t_[j] == 0.0 )
               q_[j] = 0.0;
            else if( t_[j] < 0.0 )
               q_[j] = -1.0;
            else
               q_[j] = +1.0;
      }
   }
   //
   // operator
   const Vector& operator()(const Vector& x) override
   {  //
      // n_arg, n_other
      size_t n_arg   = option_.n_arg;
      size_t n_other = option_.n_other;
      //
      // sumsq
      scalar_type sumsq(0.0);
      for(size_t j = 0; j < n_other; ++j)
      {  scalar_type model(0.0);
         scalar_type tij(1.0);
         for(size_t i = 0; i < n_arg; ++i)
         {  model += x[i] * tij;
            tij   *= t_[j];
         }
         scalar_type residual = model - q_[j];
         sumsq               += residual * residual;
      }
      //
      // y_
      y_[0] = 0.5 * sumsq;
      //
      return y_;
   }
};

} // END cmapd namespace