QuadraticCombination.cc

Go to the documentation of this file.

/*
 * $Id: QuadraticCombination.cc 177 2013-07-02 16:59:05Z jdl3 $
 * Copyright (C) 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"QuadraticCombination.hpp"

using namespace ipo_function::concrete;

QuadraticCombination::QuadraticCombination( size_t const size )
  : FunctionBase { size }, coefficients { size, size }{ coefficients.set_all( 0.0 ); }

void
QuadraticCombination::resize( size_t const size ){
  size_t const OLDSIZE { coefficients.size1() };
  if( OLDSIZE != size ){
    size_t const MIN { std::min( size, OLDSIZE ) };
    auto oldCoefficients = coefficients;
    coefficients = gsl::matrix { size, size };
    // copy old values
    for( size_t row { 0 }; row < MIN; ++row )
      for( size_t col { 0 }; col < MIN; ++col )
    coefficients.set( row, col, oldCoefficients.get( row, col ) );
    // zero new values
    for( size_t row { MIN }; row < size; ++row )
      for( size_t col { 0 }; col < size; ++col ){
    coefficients.set( row, col, 0 );
    if( row != col ) coefficients.set( col, row, 0 );
      }
  }
}

void
QuadraticCombination::setCoefficient( size_t const index, double const value ){
  if( index > coefficients.size1() ){
    using namespace ipo;
    throw IPOE( "ipo_function::concrete::QuadraticCombination::setCoefficient(): "
        "index out of range." );
  }
  coefficients.set( index, index, value );
}

void
QuadraticCombination::setCoefficient( size_t const row, size_t const col, double const value ){
  if( row == col ){
    setCoefficient( row, value );
    return;
  }
  if( row > coefficients.size1() ){
    using namespace ipo;
    throw IPOE( "ipo_function::concrete::QuadraticCombination::setCoefficient(): "
        "row out of range." );
  }
  if( col > coefficients.size2() ){
    using namespace ipo;
    throw IPOE( "ipo_function::concrete::QuadraticCombination::setCoefficient(): "
        "col out of range." );
  }
  coefficients.set( row, col, value / 2.0 );
  coefficients.set( col, row, value / 2.0 );
}

double
QuadraticCombination::getCoefficient( size_t const index ) const {
  if( index > coefficients.size1() ){
    using namespace ipo;
    throw IPOE( "ipo_function::concrete::QuadraticCombination::getCoefficient(): "
        "index out of range." );
  }
  return coefficients.get( index, index );
}

double
QuadraticCombination::getCoefficient( size_t const row, size_t const col ) const {
  if( row == col ) return getCoefficient( row );
  if( row > coefficients.size1() ){
    using namespace ipo;
    throw IPOE( "ipo_function::concrete::QuadraticCombination::getCoefficient(): "
        "row out of range." );
  }
  if( col > coefficients.size2() ){
    using namespace ipo;
    throw IPOE( "ipo_function::concrete::QuadraticCombination::getCoefficient(): "
        "col out of range." );
  }
  return 2 * coefficients.get( row, col );
}

gsl::matrix const&
QuadraticCombination::getCoefficients() const {
  return coefficients;
}

void
QuadraticCombination::setCoefficients( gsl::matrix const& matrix ){
  if( matrix.size1() != coefficients.size2() ){
    using namespace ipo;
    std::ostringstream
      ost { "ipo_function::concrete::QuadraticCombination::setCoefficients(): " };
    ost << "number of rows " << matrix.size1();
    ost << " and number of columns " << matrix.size2();
    ost << " of matrix do not match.";
    throw IPOE( ost.str() );
  }
  if( matrix.size1() != coefficients.size1() ){
    using namespace ipo;
    std::ostringstream
      ost { "ipo_function::concrete::QuadraticCombination::setCoefficients(): " };
    ost << "matrix size " << matrix.size1();
    ost << " and coefficients size " << coefficients.size1();
    ost << " do not match.";
    throw IPOE( ost.str() );
  }
  for( size_t row { 0 }; row < matrix.size1(); ++row ){
    for( size_t col { row }; col < matrix.size1(); ++col ){
      if( matrix.get( row, col ) != matrix.get( col, row ) ){
    using namespace ipo;
    std::ostringstream
      ost { "ipo_function::concrete::QuadraticCombination::setCoefficients(): " };
    ost << "matrix is not symmetric.";
    throw IPOE( ost.str() );
      }
    }
  }
  coefficients.memcpy( matrix );
}

double
QuadraticCombination::operator()( gsl::vector const& vector ){
  if( vector.size() != coefficients.size1() ){
    using namespace ipo;
    std::ostringstream ost { "ipo_function::concrete::QuadraticCombination::operator()(): " };
    ost << "vector size " << vector.size();
    ost << " and coefficients size " << coefficients.size1();
    ost << " do not match.";
    throw IPOE( ost.str() );
  }
  double result { 0.0 };
  gsl::vector temp( vector.size() );
  gsl::blas::dsymv( gsl::blas::Upper, 1, coefficients, vector, 0, temp );
  gsl::blas::ddot( temp, vector, &result );
  return result;
}

gsl::vector
QuadraticCombination::gradient( gsl::vector const& vector ){
  size_t const SIZE = vector.size();
  if( SIZE != coefficients.size1() ){
    using namespace ipo;
    std::ostringstream ost( "ipo_function::concrete::QuadraticCombination::gradient(): " );
    ost << "vector size " << vector.size();
    ost << " and coefficients size " << coefficients.size1();
    ost << " do not match.";
    throw IPOE( ost.str() );
  }
  // Computer the gradient
  if( functionGradient.size() != SIZE )
    functionGradient = gsl::vector( SIZE );
  for( size_t i { 0 }; i < SIZE; ++ i ){
    auto const row = coefficients.row( i );
    double result { 0.0 };
    gsl::blas::ddot( row, vector, &result );
    functionGradient[i] = 2.0 * result;
  }
  return functionGradient;
}

gsl::matrix
QuadraticCombination::hessian( gsl::vector const& vector ){
  size_t const SIZE = vector.size();
  if( SIZE != coefficients.size1() ){
    using namespace ipo;
    std::ostringstream ost( "ipo_function::concrete::QuadraticCombination::hessian(): " );
    ost << "vector size " << vector.size();
    ost << " and coefficients size " << coefficients.size1();
    ost << " do not match.";
    throw IPOE( ost.str() );
  }
  if( functionHessian.size1() != SIZE or functionHessian.size2() != SIZE )
    functionHessian = gsl::matrix { SIZE, SIZE };
  functionHessian.memcpy( coefficients );
  functionHessian.scale( 2.0 );
  return functionHessian;
}

void
QuadraticCombination::setVector( gsl::vector const& vector ){
  size_t const SIZE = vector.size();
  if( SIZE != coefficients.size1() ){
    using namespace ipo;
    std::ostringstream ost( "ipo_function::concrete::QuadraticaticCombination::setVector(): " );
    ost << "vector size " << vector.size();
    ost << " and coefficients size " << coefficients.size1();
    ost << " do not match.";
    throw IPOE( ost.str() );
  }
  gsl::vector temp( vector.size() );
  gsl::blas::dsymv( gsl::blas::Upper, 1.0, coefficients, vector, 0.0, temp );
  gsl::blas::ddot( temp, vector, &functionValue ); 
  if( functionGradient.size() != SIZE )
    functionGradient = gsl::vector( SIZE );
  for( size_t i { 0 }; i < SIZE; ++ i ){
    auto const row = coefficients.row( i );
    double result { 0.0 };
    gsl::blas::ddot( row, vector, &result );
    functionGradient[i] = 2.0 * result;
  }
  if( functionHessian.size1() != SIZE or functionHessian.size2() != SIZE )
    functionHessian = gsl::matrix { SIZE, SIZE };
  functionHessian.memcpy( coefficients );
  functionHessian.scale( 2.0 );
}