PhaseIFunctionAndDerivatives.cc

Go to the documentation of this file.

/*
 * $Id: PhaseIFunctionAndDerivatives.cc 236 2014-11-30 14:09:06Z jdl3 $
 * Copyright (C) 2013–2014 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"PhaseIFunctionAndDerivatives.hpp"

using namespace ipo_function::detail;

PhaseIFunctionAndDerivatives::PhaseIFunctionAndDerivatives( Function& function, double const b,
                                bool const upperBound )
  : FunctionBase { 0 == function.getSize() ? 0 : function.getSize() + 1 },
  function( function ), b { b }, upperBound { upperBound }{
  if( size != 0 ){
    functionGradient = gsl::vector( size );
    functionHessian = gsl::matrix { size, size };
  }
}

double
PhaseIFunctionAndDerivatives::operator()( gsl::vector const& vector ){
  size_t const SIZE { getSize() };
  size_t const VSIZE { vector.size() };
  
  if( 0 == VSIZE ) return 0; // do something sane
  if( 0 != SIZE and SIZE != VSIZE ) return 0; // do something sane.

  // We can now assume the first VSIZE elements of vector should
  // be passed to function
  auto subvector = vector.subvector( 0, VSIZE - 1 );
  double const& s { vector[VSIZE - 1] };
  double const f { function( subvector ) };
  if( upperBound )
    return f - b - s;
  else
    return b - f - s;
}

gsl::vector
PhaseIFunctionAndDerivatives::gradient( gsl::vector const& vector ){
  size_t const SIZE { getSize() };
  size_t const VSIZE { vector.size() };
  
  if( not functionGradient )
    functionGradient = gsl::vector( VSIZE ); // Should happen rarely

  if( 0 == VSIZE ){
    functionGradient.set_all( 0 );
    return functionGradient; // do something sane.
  }
  if( 0 != SIZE and SIZE != VSIZE ){
    functionGradient.set_all( 0 );
    return functionGradient; // do something sane.
  }
  if( 0 == SIZE and functionGradient.size() != VSIZE )
    functionGradient = gsl::vector( VSIZE ); // Should happen rarely
  
  // We can now assume the first VSIZE - 1 elements of vector should
  // be passed to gradient
  auto subvector = vector.subvector( 0, VSIZE - 1 );
  auto f = function.gradient( subvector );
  for( size_t i { 0 }; i + 1 < VSIZE; ++i )
    functionGradient[i] = upperBound ? f[i] : -f[i];
  functionGradient[VSIZE - 1] = -1;
  return functionGradient;
}

gsl::matrix
PhaseIFunctionAndDerivatives::hessian( gsl::vector const& vector ){
  size_t const SIZE { getSize() };
  size_t const VSIZE { vector.size() };
  
  if( not functionHessian )
    functionHessian = gsl::matrix { VSIZE, VSIZE }; // Should happen rarely

  if( 0 == VSIZE ){
    functionHessian.set_all( 0 );
    return functionHessian; // do something sane
  }
  if( 0 != SIZE and SIZE != VSIZE ){
    functionHessian.set_all( 0 );
    return functionHessian; // do something sane.
  }
  if( 0 == SIZE and (functionHessian.size1() != VSIZE
             or functionHessian.size2()  != VSIZE) )
    functionHessian = gsl::matrix { VSIZE, VSIZE }; // Should happen rarely
  
  // We can now assume the first VSIZE - 1 elements of vector should
  // be passed to hessian
  auto subvector = vector.subvector( 0, VSIZE - 1 );
  auto h = function.hessian( subvector );
  if( upperBound ){
    for( size_t i { 0 }; i < VSIZE; ++i )
      for( size_t j { i }; j < VSIZE; ++j ){
    double const value { (i + 1 == VSIZE or j + 1 == VSIZE) ? 0 : h.get( i, j ) };
    functionHessian.set( i, j, value );
    if( j > i ) functionHessian.set( j, i, value );
      }
  } else {
    for( size_t i { 0 }; i < VSIZE; ++i )
      for( size_t j { i }; j < VSIZE; ++j ){
    double const value { (i + 1 == VSIZE or j + 1 == VSIZE) ? 0 : -h.get( i, j ) };
    functionHessian.set( i, j, value );
    if( j > i ) functionHessian.set( j, i, value );
      }
  }
  return functionHessian;
}

void
PhaseIFunctionAndDerivatives::setVector( gsl::vector const& vector ){
  size_t const SIZE { getSize() };
  size_t const VSIZE { vector.size() };

  if( not functionGradient )
    functionGradient = gsl::vector( VSIZE ); // Should happen rarely
  if( not functionHessian )
    functionHessian = gsl::matrix( VSIZE, VSIZE ); // Should happen rarely

  if( 0 == VSIZE ){
    functionGradient.set_all( 0 );
    functionHessian.set_all( 0 );
    return; // do something sane
  }
  if( 0 != SIZE and SIZE != VSIZE ){
    functionGradient.set_all( 0 );
    functionHessian.set_all( 0 );
    return; // do something sane
  }
  if( 0 == SIZE and functionGradient.size() != VSIZE )
    functionGradient = gsl::vector( VSIZE ); // Should happen rarely
  if( 0 == SIZE and (functionHessian.size1() != VSIZE
             or functionHessian.size2() != VSIZE) )
    functionHessian = gsl::matrix { VSIZE, VSIZE }; // Should happen rarely
  
  // Deal separately with cases where function is or is not a DerivativesEstimates object
  auto de = dynamic_cast<ipo_function::DerivativesEstimates*>( &function );
  if( nullptr == de ){ // We need to estimate gradient and hessian as best we can
    functionValue = (*this)( vector );
    gradient( vector );
    hessian( vector );
  } else { // We are dealing with a DerivativesEstimates object
    auto subvector = vector.subvector( 0, VSIZE - 1 );
    de->setVector( subvector );
    double const& s { vector[VSIZE - 1] };
    if( upperBound )
      functionValue = de->value() - b - s;
    else
      functionValue = b - de->value() - s;
    for( size_t i { 0 }; i < VSIZE; ++i ){
      auto fG = de->gradient();
      auto fH = de->hessian();
      if( i + 1 < VSIZE ) functionGradient[i] = upperBound ? fG[i] : -fG[i];
      else functionGradient[i] = -1;
      if( upperBound ){
    for( size_t j { i }; j < VSIZE; ++j ){
      double const h { (i + 1 == VSIZE or j + 1 == VSIZE) ? 0 : fH.get( i, j ) };
      functionHessian.set( i, j, h );
      if( j > i ) functionHessian.set( j, i, h );
    }
      } else {
    for( size_t j { i }; j < VSIZE; ++j ){
      double const h { (i + 1 == VSIZE or j + 1 == VSIZE) ? 0 : -fH.get( i, j ) };
      functionHessian.set( i, j, h );
      if( j > i ) functionHessian.set( j, i, h );
    }
      }
    }
  }
  
}