Function.cc

Go to the documentation of this file.

/*
 * $Id: Function.cc 139 2013-06-29 15:11:33Z jdl3 $
 * Copyright (C) 2011, 2012, 2013 John D Lamb
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#ifdef HAVE_CONFIG_H
#  include<config.h>
#endif

#include"Function.hpp"

using namespace ipo_function;

Function::Function( size_t const size ) : FunctionBase { size }{}

gsl::vector
Function::gradient( gsl::vector const& vector ){
  // This is where all the work is done. 
  double const H { 2 * h };
  // Find SIZE
  const size_t SIZE { vector.size() };
  // Set up gradient
  auto gradient = gsl::vector( SIZE );
  // Fill diffVector
  for( size_t i { 0 }; i < SIZE; ++i ){
    double tmp { vector.get( i ) };
    const_cast<gsl::vector&>( vector ).set( i, tmp + h );
    double const u { operator()( vector ) };
    const_cast<gsl::vector&>( vector ).set( i, tmp - h );
    double const l { operator()( vector ) };
    const_cast<gsl::vector&>( vector ).set( i, tmp );
    //std::cout << "(" << l << "," << u << "): " << ((u - l) / H) << std::endl;
    // vector returned to original value.
    gradient.set( i, (u - l) / H );
  }
  return gradient;
}

gsl::matrix
Function::hessian( gsl::vector const& vector ){
  // This is where all the work is done. 
  double const H { 4 * h2 * h2 };
  //std::cout << "h = " << h << std::endl;
  //std::cout << "H = " << H << std::endl;
  // Find SIZE
  const size_t SIZE { vector.size() };
  // Set up Hessian
  auto hessian = gsl::matrix( SIZE, SIZE );
  for( size_t i { 0 }; i < SIZE; ++i ){
    // diagonal elements
    double const c { operator()( vector ) };
    double tmpi { vector.get( i ) };
    const_cast<gsl::vector&>( vector ).set( i, tmpi + 2 * h2 );
    double const a { operator()( vector ) };
    const_cast<gsl::vector&>( vector ).set( i, tmpi - 2 * h2 );
    double const d { operator()( vector ) };
    const_cast<gsl::vector&>( vector ).set( i, tmpi );
    hessian.set( i, i, (a - 2 * c + d) / H );
    for( size_t j { i + 1 }; j < SIZE; ++j ){
      // off diagonal elements
      double tmpj { vector.get( j ) };
      const_cast<gsl::vector&>( vector ).set( i, tmpi + h2 );
      const_cast<gsl::vector&>( vector ).set( j, tmpj + h2 );
      double const a { operator()( vector ) };
      const_cast<gsl::vector&>( vector ).set( j, tmpj - h2 );
      double const b { operator()( vector ) };
      const_cast<gsl::vector&>( vector ).set( i, tmpi - h2 );
      double const d { operator()( vector ) };
      const_cast<gsl::vector&>( vector ).set( j, tmpj + h2 );
      double const c { operator()( vector ) };
      const_cast<gsl::vector&>( vector ).set( i, tmpi );
      const_cast<gsl::vector&>( vector ).set( j, tmpj );
      //std::cout << "a = " << a << std::endl;
      //std::cout << "b = " << b << std::endl;
      //std::cout << "c = " << c << std::endl;
      //std::cout << "d = " << d << std::endl;
      double const entry = (a - b - c + d) / H;
      hessian.set( i, j, entry );
      hessian.set( j, i, entry );
    }
  }
  // DEBUG
  // std::cout << "Hessian (estimated by central differences)" << std::endl;
  // for( auto row : hessian ){
  //   for( auto entry : row )
  //     std::cout << entry << " ";
  //   std::cout << std::endl;
  // }
  // std::cout << "-----------------------" << std::endl;
  // END DEBUG
  return hessian;
}

std::tuple<gsl::vector,gsl::matrix>
Function::derivatives( gsl::vector const& vector ){
  return std::make_tuple( gradient( vector ), hessian( vector ) );
}