problem/base.cpp

/*****************************************************************************
 *   Copyright (C) 2004-2015 The PaGMO development team,                     *
 *   Advanced Concepts Team (ACT), European Space Agency (ESA)               *
 *                                                                           *
 *   https://github.com/esa/pagmo                                            *
 *                                                                           *
 *   act@esa.int                                                             *
 *                                                                           *
 *   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.,                                         *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.               *
 *****************************************************************************/

// 30/01/10 Created by Francesco Biscani.

#include <algorithm>
#include <boost/math/special_functions/fpclassify.hpp>
#include <boost/numeric/conversion/bounds.hpp>
#include <boost/numeric/conversion/cast.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/ref.hpp>
#include <cmath>
#include <climits>
#include <cstddef>
#include <iostream>
#include <iterator>
#include <numeric>
#include <sstream>
#include <string>
#include <typeinfo>
#include <vector>

#include "../exceptions.h"
#include "../population.h"
#include "../types.h"
#include "base.h"

namespace pagmo
{
namespace problem {


base::base(int n, int ni, int nf, int nc, int nic, const double &c_tol): //TODO should we use size_type directly?
        m_i_dimension(boost::numeric_cast<size_type>(ni)),m_f_dimension(boost::numeric_cast<f_size_type>(nf)),
        m_c_dimension(boost::numeric_cast<c_size_type>(nc)),m_ic_dimension(boost::numeric_cast<c_size_type>(nic)),
        m_c_tol(nc,c_tol),
        m_decision_vector_cache_f(boost::numeric_cast<decision_vector_cache_type::size_type>(cache_capacity)),
        m_fitness_vector_cache(boost::numeric_cast<fitness_vector_cache_type::size_type>(cache_capacity)),
        m_decision_vector_cache_c(boost::numeric_cast<decision_vector_cache_type::size_type>(cache_capacity)),
        m_constraint_vector_cache(boost::numeric_cast<constraint_vector_cache_type::size_type>(cache_capacity)),
        m_best_x(0),
        m_best_f(0),
        m_best_c(0),
        m_fevals(0),
        m_cevals(0)
{
        if (c_tol < 0) {
                pagmo_throw(value_error,"constraints tolerance must be non-negative");
        }
        if (n <= 0 || !nf || ni > n || nic > nc) {
                pagmo_throw(value_error,"invalid dimension(s)");
        }
        const size_type size = boost::numeric_cast<size_type>(n);
        m_lb.resize(size);
        m_ub.resize(size);
        std::fill(m_lb.begin(),m_lb.end(),0);
        std::fill(m_ub.begin(),m_ub.end(),1);
        // Resize properly temporary fitness and constraint storage.
        m_tmp_f1.resize(m_f_dimension);
        m_tmp_f2.resize(m_f_dimension);
        m_tmp_c1.resize(m_c_dimension);
        m_tmp_c2.resize(m_c_dimension);
        // Normalise bounds.
        normalise_bounds();
}


base::base(int n, int ni, int nf, int nc, int nic, const std::vector<double> &c_tol): //TODO should we use size_type directly?
        m_i_dimension(boost::numeric_cast<size_type>(ni)),m_f_dimension(boost::numeric_cast<f_size_type>(nf)),
        m_c_dimension(boost::numeric_cast<c_size_type>(nc)),m_ic_dimension(boost::numeric_cast<c_size_type>(nic)),
        m_c_tol(c_tol),
        m_decision_vector_cache_f(boost::numeric_cast<decision_vector_cache_type::size_type>(cache_capacity)),
        m_fitness_vector_cache(boost::numeric_cast<fitness_vector_cache_type::size_type>(cache_capacity)),
        m_decision_vector_cache_c(boost::numeric_cast<decision_vector_cache_type::size_type>(cache_capacity)),
        m_constraint_vector_cache(boost::numeric_cast<constraint_vector_cache_type::size_type>(cache_capacity)),
        m_best_x(0),
        m_best_f(0),
        m_best_c(0),
        m_fevals(0),
        m_cevals(0)
{
        if (c_tol.size() != static_cast<constraint_vector::size_type>(nc) ) {
                pagmo_throw(value_error,"invalid constraints vector dimension");
        }
        for(unsigned int i=0; i<(m_c_dimension - m_ic_dimension);i++) {
                if (c_tol[i] < 0) {
                        pagmo_throw(value_error,"constraints tolerance must be non-negative for equality constraints");
                }
        }
        if (n <= 0 || !nf || ni > n || nic > nc) {
                pagmo_throw(value_error,"invalid dimension(s)");
        }
        const size_type size = boost::numeric_cast<size_type>(n);
        m_lb.resize(size);
        m_ub.resize(size);
        std::fill(m_lb.begin(),m_lb.end(),0);
        std::fill(m_ub.begin(),m_ub.end(),1);
        // Resize properly temporary fitness and constraint storage.
        m_tmp_f1.resize(m_f_dimension);
        m_tmp_f2.resize(m_f_dimension);
        m_tmp_c1.resize(m_c_dimension);
        m_tmp_c2.resize(m_c_dimension);
        // Normalise bounds.
        normalise_bounds();
}


base::base(const double &l_value, const double &u_value, int n, int ni, int nf, int nc, int nic, const double &c_tol):
        m_i_dimension(boost::numeric_cast<size_type>(ni)),m_f_dimension(boost::numeric_cast<f_size_type>(nf)),
        m_c_dimension(boost::numeric_cast<c_size_type>(nc)),m_ic_dimension(boost::numeric_cast<c_size_type>(nic)),
        m_c_tol(nc,c_tol),
        m_decision_vector_cache_f(boost::numeric_cast<decision_vector_cache_type::size_type>(cache_capacity)),
        m_fitness_vector_cache(boost::numeric_cast<fitness_vector_cache_type::size_type>(cache_capacity)),
        m_decision_vector_cache_c(boost::numeric_cast<decision_vector_cache_type::size_type>(cache_capacity)),
        m_constraint_vector_cache(boost::numeric_cast<constraint_vector_cache_type::size_type>(cache_capacity)),
        m_best_x(0),
        m_best_f(0),
        m_best_c(0),
        m_fevals(0),
        m_cevals(0)
{
        if (c_tol < 0) {
                pagmo_throw(value_error,"constraints tolerance must be non-negative");
        }
        if (n <= 0 || !nf || ni > n || nic > nc) {
                pagmo_throw(value_error,"invalid dimension(s)");
        }
        if (l_value > u_value) {
                pagmo_throw(value_error,"value for lower bounds cannot be greater than value for upper bounds");
        }
        const size_type size = boost::numeric_cast<size_type>(n);
        m_lb.resize(size);
        m_ub.resize(size);
        std::fill(m_lb.begin(),m_lb.end(),l_value);
        std::fill(m_ub.begin(),m_ub.end(),u_value);
        // Resize properly temporary fitness and constraint storage.
        m_tmp_f1.resize(m_f_dimension);
        m_tmp_f2.resize(m_f_dimension);
        m_tmp_c1.resize(m_c_dimension);
        m_tmp_c2.resize(m_c_dimension);
        // Normalise bounds.
        normalise_bounds();
}


base::base(const decision_vector &lb, const decision_vector &ub, int ni, int nf, int nc, int nic, const double &c_tol):
        m_i_dimension(boost::numeric_cast<size_type>(ni)),m_f_dimension(boost::numeric_cast<f_size_type>(nf)),
        m_c_dimension(boost::numeric_cast<c_size_type>(nc)),m_ic_dimension(boost::numeric_cast<c_size_type>(nic)),
        m_c_tol(nc,c_tol),
        m_decision_vector_cache_f(boost::numeric_cast<decision_vector_cache_type::size_type>(cache_capacity)),
        m_fitness_vector_cache(boost::numeric_cast<fitness_vector_cache_type::size_type>(cache_capacity)),
        m_decision_vector_cache_c(boost::numeric_cast<decision_vector_cache_type::size_type>(cache_capacity)),
        m_constraint_vector_cache(boost::numeric_cast<constraint_vector_cache_type::size_type>(cache_capacity)),
        m_best_x(0),
        m_best_f(0),
        m_best_c(0),
        m_fevals(0),
        m_cevals(0)
{
        if (c_tol < 0) {
                pagmo_throw(value_error,"constraints tolerance must be non-negative");
        }
        if (!nf || m_i_dimension > lb.size() || nic > nc) {
                pagmo_throw(value_error,"invalid dimension(s)");
        }
        construct_from_iterators(lb.begin(),lb.end(),ub.begin(),ub.end());
        // Resize properly temporary fitness and constraint storage.
        m_tmp_f1.resize(m_f_dimension);
        m_tmp_f2.resize(m_f_dimension);
        m_tmp_c1.resize(m_c_dimension);
        m_tmp_c2.resize(m_c_dimension);
        // Normalise bounds.
        normalise_bounds();
}


base::~base() {}


std::string base::get_name() const
{
        return typeid(*this).name();
}


const decision_vector &base::get_lb() const
{
        return m_lb;
}


const decision_vector &base::get_ub() const
{
        return m_ub;
}


void base::set_bounds(const decision_vector &lb, const decision_vector &ub)
{
        if (lb.size() != ub.size() || lb.size() != m_lb.size()) {
                pagmo_throw(value_error,"invalid or inconsistent bounds dimensions in set_bounds()");
        }
        verify_bounds(lb.begin(),lb.end(),ub.begin(),ub.end());
        m_lb = lb;
        m_ub = ub;
        // Normalise bounds.
        normalise_bounds();
}


void base::set_bounds(const double &l_value, const double &u_value)
{
        if (l_value > u_value) {
                pagmo_throw(value_error,"lower bound cannot be greater than upper bound in set_bounds()");
        }
        std::fill(m_lb.begin(),m_lb.end(),l_value);
        std::fill(m_ub.begin(),m_ub.end(),u_value);
        normalise_bounds();
}


void base::set_bounds(int n, const double &l_value, const double &u_value)
{
        if (l_value > u_value) {
                pagmo_throw(value_error,"lower bound cannot be greater than upper bound in set_bounds()");
        }
        m_lb[n] = l_value;
        m_ub[n] = u_value;
        // Normalise bounds.
        normalise_bounds();
}


void base::set_lb(const decision_vector &lb)
{
        if (lb.size() != m_lb.size()) {
                pagmo_throw(value_error,"invalid bounds dimension in set_lb()");
        }
        verify_bounds(lb.begin(),lb.end(),m_ub.begin(),m_ub.end());
        m_lb = lb;
        // Normalise bounds.
        normalise_bounds();
}


void base::set_lb(int n, const double &value)
{
        const size_type i = boost::numeric_cast<size_type>(n);
        if (i >= m_lb.size() || m_ub[i] < value) {
                pagmo_throw(value_error,"invalid index and/or value for lower bound");
        }
        m_lb[i] = value;
        // Normalise bounds.
        normalise_bounds();
}


void base::set_lb(const double &value)
{
        for (size_type i = 0; i < m_lb.size(); ++i) {
                if (m_ub[i] < value) {
                        pagmo_throw(value_error,"invalid value for lower bound");
                }
        }
        std::fill(m_lb.begin(),m_lb.end(),value);
        // Normalise bounds.
        normalise_bounds();
}


void base::set_ub(const decision_vector &ub)
{
        if (ub.size() != m_lb.size()) {
                pagmo_throw(value_error,"invalid bounds dimension in set_ub()");
        }
        verify_bounds(m_lb.begin(),m_lb.end(),ub.begin(),ub.end());
        m_ub = ub;
        // Normalise bounds.
        normalise_bounds();
}


void base::set_ub(int n, const double &value)
{
        const size_type i = boost::numeric_cast<size_type>(n);
        if (i >= m_lb.size() || m_lb[i] > value) {
                pagmo_throw(value_error,"invalid index and/or value for upper bound");
        }
        m_ub[i] = value;
        // Normalise bounds.
        normalise_bounds();
}


void base::set_ub(const double &value)
{
        for (size_type i = 0; i < m_lb.size(); ++i) {
                if (m_lb[i] > value) {
                        pagmo_throw(value_error,"invalid value for upper bound");
                }
        }
        std::fill(m_ub.begin(),m_ub.end(),value);
        // Normalise bounds.
        normalise_bounds();
}


unsigned int base::get_fevals() const
{
        return m_fevals;
}


unsigned int base::get_cevals() const
{
        return m_cevals;
}



base::size_type base::get_dimension() const
{
        return m_lb.size();
}


base::size_type base::get_i_dimension() const
{
        return m_i_dimension;
}


base::f_size_type base::get_f_dimension() const
{
        return m_f_dimension;
}


base::c_size_type base::get_c_dimension() const
{
        return m_c_dimension;
}


base::c_size_type base::get_ic_dimension() const
{
        return m_ic_dimension;
}


const std::vector<double>& base::get_c_tol() const
{
        return m_c_tol;
}


double base::get_diameter() const
{
        double retval = 0;
        for (size_type i = 0; i < get_dimension(); ++i) {
                retval += (m_ub[i] - m_lb[i]) * (m_ub[i] - m_lb[i]);
        }
        return std::sqrt(retval);
}


fitness_vector base::objfun(const decision_vector &x) const
{
        if (x.size()!=get_dimension()) {
                pagmo_throw(value_error,"invalid chromosome length");
        }
        fitness_vector f(m_f_dimension);
        objfun(f,x);
        return f;
}


void base::objfun(fitness_vector &f, const decision_vector &x) const
{
        // Some checks on the input values.
        if (f.size() != m_f_dimension) {
                pagmo_throw(value_error,"wrong fitness vector size when calling objective function");
        }
        if (x.size() != get_dimension()) {
                pagmo_throw(value_error,"wrong decision vector size when calling objective function");
        }
        // Look into the cache.
        typedef decision_vector_cache_type::iterator x_iterator;
        typedef fitness_vector_cache_type::iterator f_iterator;
        const x_iterator x_it = std::find(m_decision_vector_cache_f.begin(),m_decision_vector_cache_f.end(),x);
        if (x_it == m_decision_vector_cache_f.end()) {
                // Fitness is not into memory. Calculate it.
                objfun_impl(f,x);
                // Increase function evaluation counter.
                m_fevals++;
                // Make sure that the implementation of objfun_impl() in the derived class did not fuck up the dimension of the fitness vector.
                if (f.size() != m_f_dimension) {
                        pagmo_throw(value_error,"fitness dimension was changed inside objfun_impl()");
                }
                // Store the decision vector and the newly-calculated fitness in the front of the buffers.
                m_decision_vector_cache_f.push_front(x);
                m_fitness_vector_cache.push_front(f);
        } else {
                // Compute the corresponding iterator in the fitness vector cache.
                f_iterator f_it = m_fitness_vector_cache.begin();
                std::advance(f_it,std::distance(m_decision_vector_cache_f.begin(),x_it));
                pagmo_assert(f_it != m_fitness_vector_cache.end());
                // Assign to the fitness vector the value in the cache.
                f = *f_it;
                // Move the content of the current positions to the front of the buffers, and shift everything else down
                // by one position.
                x_iterator tmp_x_it = m_decision_vector_cache_f.begin();
                f_iterator tmp_f_it = m_fitness_vector_cache.begin();
                while (x_it != tmp_x_it) {
                        x_it->swap(*tmp_x_it);
                        f_it->swap(*tmp_f_it);
                        ++tmp_x_it;
                        ++tmp_f_it;
                }
                pagmo_assert(tmp_f_it == f_it);
        }
}


bool base::compare_fitness(const fitness_vector &v_f1, const fitness_vector &v_f2) const
{
        if (v_f1.size() != m_f_dimension || v_f2.size() != m_f_dimension) {
                pagmo_throw(value_error,"invalid sizes for fitness vector(s) during comparison");
        }
        return compare_fitness_impl(v_f1,v_f2);
}


bool base::compare_fitness_impl(const fitness_vector &v_f1, const fitness_vector &v_f2) const
{
        pagmo_assert(v_f1.size() == v_f2.size()); // NOTE: v_f1.size CAN be different from m_f_dimension (i.e. when its a base_meta calling it)
        f_size_type count1 = 0; f_size_type count2 = 0;
        for (f_size_type i = 0; i < v_f1.size(); ++i) {
                if (v_f1[i] < v_f2[i]) {
                        ++count1;
                }
                if (v_f1[i] == v_f2[i]) {
                        ++count2;
                }
        }
        return ( ( (count1+count2) == v_f1.size()) && (count1>0) );
}


std::string base::human_readable() const
{
        std::ostringstream s;
        s << "Problem name: " << get_name() << '\n';
        const size_type size = get_dimension();
        s << "\tGlobal dimension:\t\t\t" << size << '\n';
        s << "\tInteger dimension:\t\t\t" << m_i_dimension << '\n';
        s << "\tFitness dimension:\t\t\t" << m_f_dimension << '\n';
        s << "\tConstraints dimension:\t\t\t" << m_c_dimension << '\n';
        s << "\tInequality constraints dimension:\t" << m_ic_dimension << '\n';
        s << "\tLower bounds: ";
        s << m_lb << '\n';
        s << "\tUpper bounds: ";
        s << m_ub << '\n';
        s << "\tConstraints tolerance: ";
        s << m_c_tol << '\n';
        s << human_readable_extra();
        return s.str();
}


std::string base::human_readable_extra() const
{
        return std::string();
}


bool base::operator==(const base &p) const
{
        if (!is_compatible(p)) {
                return false;
        }
        return equality_operator_extra(p);
}


bool base::is_compatible(const base &p) const
{
// marcusm: changed this for experiments with migration between single- and multiobjective islands.
/*      if (typeid(*this) != typeid(p) || get_dimension() != p.get_dimension() || m_i_dimension != p.m_i_dimension || m_f_dimension != p.m_f_dimension ||
                m_c_dimension != p.m_c_dimension || m_ic_dimension != p.m_ic_dimension)*/
        if (typeid(*this) != typeid(p) || get_dimension() != p.get_dimension() || m_i_dimension != p.m_i_dimension ||
                m_c_dimension != p.m_c_dimension || m_ic_dimension != p.m_ic_dimension)
        {
                return false;
        }
        return true;
}


bool base::compare_x(const decision_vector &x1, const decision_vector &x2) const
{
        // Make sure the size of the tmp fitness vectors are suitable.
        pagmo_assert(m_tmp_f1.size() == m_f_dimension && m_tmp_f2.size() == m_f_dimension);
        // Store fitnesses into temporary space.
        objfun(m_tmp_f1,x1);
        objfun(m_tmp_f2,x2);
        // Make sure the size of the tmp constraint vectors are suitable.
        pagmo_assert(m_tmp_c1.size() == m_c_dimension && m_tmp_c2.size() == m_c_dimension);
        // Store constraint vector into temporary space.
        compute_constraints(m_tmp_c1,x1);
        compute_constraints(m_tmp_c2,x2);
        // Call the comparison implementation.
        return compare_fc(m_tmp_f1,m_tmp_c1,m_tmp_f2,m_tmp_c2);
}


bool base::compare_fc(const fitness_vector &f1, const constraint_vector &c1, const fitness_vector &f2, const constraint_vector &c2) const
{
        if (f1.size() != m_f_dimension || f2.size() != m_f_dimension) {
                pagmo_throw(value_error,"wrong size(s) for fitness vector(s)");
        }
        if (c1.size() != m_c_dimension || c2.size() != m_c_dimension) {
                pagmo_throw(value_error,"wrong size(s) for constraint vector(s)");
        }
        if (m_c_dimension) {
                return compare_fc_impl(f1,c1,f2,c2);
        } else {
                return compare_fitness_impl(f1,f2);
        }
}


bool base::compare_fc_impl(const fitness_vector &f1, const constraint_vector &c1, const fitness_vector &f2, const constraint_vector &c2) const
{
        const bool test1 = feasibility_c(c1), test2 = feasibility_c(c2);
        if (test1 && !test2) {
                return true;
        }
        if (!test1 && test2) {
                return false;
        }
        // At this point, either they are both feasible or they are not.
        if (test1) {
                pagmo_assert(test2);
                // If they are feasible, compare fitnesses and return.
                return compare_fitness_impl(f1,f2);
        } else {
                pagmo_assert(!test2);
                // If they are not feasible, compare constraints and return.
                return compare_constraints_impl(c1,c2);
        }
}



void base::compute_constraints_impl(constraint_vector &c, const decision_vector &x) const
{
        (void)x;
        std::fill(c.begin(),c.end(),0);
}


void base::compute_constraints(constraint_vector &c, const decision_vector &x) const
{
        if (x.size() != get_dimension() || c.size() != get_c_dimension()) {
                pagmo_throw(value_error,"invalid constraint and/or decision vector(s) size(s) during constraint testing");
        }
        // Do not do anything if constraints size is 0.
        if (!m_c_dimension) {
                return;
        }
        // Look into the cache.
        typedef decision_vector_cache_type::iterator x_iterator;
        typedef constraint_vector_cache_type::iterator c_iterator;
        const x_iterator x_it = std::find(m_decision_vector_cache_c.begin(),m_decision_vector_cache_c.end(),x);
        if (x_it == m_decision_vector_cache_c.end()) {
                // Constraint vector is not into memory. Calculate it.
                compute_constraints_impl(c,x);
                m_cevals++;
                // Make sure c was not fucked up in the implementation of constraints calculation.
                if (c.size() != get_c_dimension()) {
                        pagmo_throw(value_error,"constraints dimension was changed inside compute_constraints_impl()");
                }
                // Store the decision vector and the newly-calculated constraint vector in the front of the buffers.
                m_decision_vector_cache_c.push_front(x);
                m_constraint_vector_cache.push_front(c);
        } else {
                // Compute the corresponding iterator in the constraint vector cache.
                c_iterator c_it = m_constraint_vector_cache.begin();
                std::advance(c_it,std::distance(m_decision_vector_cache_c.begin(),x_it));
                pagmo_assert(c_it != m_constraint_vector_cache.end());
                // Assign to the fitness vector the value in the cache.
                c = *c_it;
                // Move the content of the current positions to the front of the buffers, and shift everything else down
                // by one position.
                x_iterator tmp_x_it = m_decision_vector_cache_c.begin();
                c_iterator tmp_c_it = m_constraint_vector_cache.begin();
                while (x_it != tmp_x_it) {
                        x_it->swap(*tmp_x_it);
                        c_it->swap(*tmp_c_it);
                        ++tmp_x_it;
                        ++tmp_c_it;
                }
                pagmo_assert(tmp_c_it == c_it);
        }
}


constraint_vector base::compute_constraints(const decision_vector &x) const
{
        constraint_vector c(get_c_dimension());
        compute_constraints(c,x);
        return c;
}


bool base::feasibility_x(const decision_vector &x) const
{
        // Compute the constraints and store internally.
        compute_constraints(m_tmp_c1,x);
        return feasibility_c(m_tmp_c1);
}


bool base::test_constraint(const constraint_vector &c, const c_size_type &i) const
{
        pagmo_assert(i < m_c_dimension);
        if (i < m_c_dimension - m_ic_dimension) {
                // Equality constraint testing.
                return (std::abs(c[i]) <= m_c_tol[i]);
        } else {
                return c[i] <= m_c_tol[i];
        }
}


bool base::feasibility_c(const constraint_vector &c) const
{
        if (c.size() != m_c_dimension) {
                pagmo_throw(value_error,"invalid size for constraint vector");
        }
        pagmo_assert(m_c_dimension >= m_ic_dimension);
        for (c_size_type i = 0; i < m_c_dimension; ++i) {
                if (!test_constraint(c,i)) {
                        return false;
                }
        }
        return true;
}


bool base::compare_constraints(const constraint_vector &c1, const constraint_vector &c2) const
{
        if (c1.size() != m_c_dimension || c2.size() != m_c_dimension) {
                pagmo_throw(value_error,"invalid size(s) for constraint vector(s)");
        }
        return compare_constraints_impl(c1,c2);
}


bool base::compare_constraints_impl(const constraint_vector &c1, const constraint_vector &c2) const
{
        pagmo_assert(c1.size() == c2.size() && c1.size() == m_c_dimension);
        // Counters of satisfied constraints.
        c_size_type count1 = 0, count2 = 0;
        // L2 norm of constraints mismatches.
        double norm1 = 0, norm2 = 0;
        // Equality constraints.
        for (c_size_type i = 0; i < m_c_dimension - m_ic_dimension; ++i) {
                if (test_constraint(c1,i)) {
                        ++count1;
                }
                if (test_constraint(c2,i)) {
                        ++count2;
                }
                norm1 += std::abs(c1[i]) * std::abs(c1[i]);
                norm2 += std::abs(c2[i]) * std::abs(c2[i]);
        }
        // Inequality constraints.
        for (c_size_type i = m_c_dimension - m_ic_dimension; i < m_c_dimension; ++i) {
                if (test_constraint(c1,i)) {
                        ++count1;
                } else {
                        norm1 += c1[i] * c1[i];
                }
                if (test_constraint(c2,i)) {
                        ++count2;
                } else {
                        norm2 += c2[i] * c2[i];
                }
        }
        if (count1 > count2) {
                return true;
        } else if (count1 < count2) {
                return false;
        } else {
                return (norm1 < norm2);
        }
}



bool base::equality_operator_extra(const base &p) const
{
        (void)p;
        return true;
}


bool base::operator!=(const base &p) const
{
        return !(*this == p);
}


bool base::verify_x(const decision_vector &x) const
{
        if (x.size() != get_dimension()) {
                return false;
        }
        for (size_type i = 0; i < get_dimension(); ++i) {
                if (x[i] < m_lb[i] || x[i] > m_ub[i]) {
                        return false;
                }
                // Check that, if this is an integer component, it is really an integer.
                if (i >= (get_dimension() - m_i_dimension) && double_to_int::convert(x[i]) != x[i]) {
                        return false;
                }
        }
        return true;
}

// This function will round to the nearest integer the upper/lower bounds of the integer part of the problem.
// This should be called each time bounds are set.
void base::normalise_bounds()
{
        pagmo_assert(m_lb.size() >= m_i_dimension);
        // Flag to be set if we had to fix the bounds.
        bool bounds_fixed = false;
        for (size_type i = 0; i < m_lb.size() - m_i_dimension; ++i) {
                // Handle NaNs: if either lower/upper bound(s) are NaN, replace with 0 and 1 respectively.
                if (boost::math::isnan(m_lb[i]) || boost::math::isnan(m_ub[i])) {
                        m_lb[i] = 0;
                        m_ub[i] = 1;
                        bounds_fixed = true;
                }
                // +-Infs are replaced by the highest/lowest values representable by double.
                if (boost::math::isinf(m_lb[i])) {
                        m_lb[i] = (m_lb[i] > 0) ? boost::numeric::bounds<double>::highest() : boost::numeric::bounds<double>::lowest();
                        bounds_fixed = true;
                }
                if (boost::math::isinf(m_ub[i])) {
                        m_ub[i] = (m_ub[i] > 0) ? boost::numeric::bounds<double>::highest() : boost::numeric::bounds<double>::lowest();
                        bounds_fixed = true;
                }
        }
        for (size_type i = m_lb.size() - m_i_dimension; i < m_lb.size(); ++i) {
                // First let's make sure that integer bounds are in the allowed range.
                if (m_lb[i] < INT_MIN) {
                        m_lb[i] = INT_MIN;
                        bounds_fixed = true;
                }
                if (m_lb[i] > INT_MAX) {
                        m_lb[i] = INT_MAX;
                        bounds_fixed = true;
                }
                if (m_ub[i] < INT_MIN) {
                        m_ub[i] = INT_MIN;
                        bounds_fixed = true;
                }
                if (m_ub[i] > INT_MAX) {
                        m_ub[i] = INT_MAX;
                        bounds_fixed = true;
                }
                // Then convert them to the nearest integer if necessary.
                if (m_lb[i] != double_to_int::convert(m_lb[i])) {
                        m_lb[i] = double_to_int::convert(m_lb[i]);
                        bounds_fixed = true;
                }
                if (m_ub[i] != double_to_int::convert(m_ub[i])) {
                        m_ub[i] = double_to_int::convert(m_ub[i]);
                        bounds_fixed = true;
                }
        }
        if (bounds_fixed) {
                pagmo_throw(value_error,"problem bounds were invalid and had to be fixed");
        }
}


std::ostream &operator<<(std::ostream &s, const base &p)
{
        s << p.human_readable();
        return s;
}


void base::set_sparsity(int &lenG, std::vector<int> &iGfun, std::vector<int> &jGvar) const
{
        (void)lenG;
        (void)iGfun;
        (void)jGvar;
        pagmo_throw(not_implemented_error,"sparsity is not implemented for this problem");
}


void base::estimate_sparsity(const decision_vector &x0, int& lenG, std::vector<int>& iGfun, std::vector<int>& jGvar) const {
        // We check that the user is providing a decision vector that is of the required length
        if (!verify_x(x0)) {
                pagmo_throw(value_error,"cannot estimate pattern from this decision vector: not compatible with problem");
        }
        size_type Dc = m_lb.size() - m_i_dimension;
        fitness_vector f0(m_f_dimension),f_new(m_f_dimension);
        objfun(f0,x0);
        constraint_vector c0(m_c_dimension),c_new(m_c_dimension);
        compute_constraints(c0,x0);
        decision_vector x_new = x0;
        iGfun.resize(0);jGvar.resize(0); lenG=0;

        for (size_type j=0;j<Dc;++j)
        {
                //we perturb the component of x0 only if ub>lb, if ub=lb the variable is assumed
                //to be 'just' a parameter ... in some problem implementations this is rather
                //useful, but it also requires that the algorithm treat those variables accordingly (i.e.
                //it does not allow them to be outside the box bounds)
                if (m_ub[j] == m_lb[j]) continue;
                x_new[j] = x0[j] +  std::max(std::fabs(x0[j]), 1.0) * 1e-8;
                objfun(f_new,x_new);
                compute_constraints(c_new,x_new);
                for (size_type i=0;i<m_f_dimension;++i)
                {
                        if (f_new[i]!=f0[i]) {iGfun.push_back(i); jGvar.push_back(j); lenG++;}
                }
                for (size_type i=0;i<m_c_dimension;++i)
                {
                        if (c_new[i]!=c0[i]) {iGfun.push_back(i+m_f_dimension); jGvar.push_back(j); lenG++;}
                }
                x_new[j] = x0[j];
        }

}


void base::estimate_sparsity(int& lenG, std::vector<int>& iGfun, std::vector<int>& jGvar) const {

        size_type Dc = m_lb.size() - m_i_dimension;
        fitness_vector f0(m_f_dimension),f_new(m_f_dimension);
        decision_vector x0(Dc);
        // Double precision random number generator.
        rng_double drng(rng_generator::get<rng_double>());

        for (decision_vector::size_type i = 0; i<Dc;++i) {
                x0[i] = boost::uniform_real<double>(m_lb[i],m_ub[i])(drng);
        }

        objfun(f0,x0);
        constraint_vector c0(m_c_dimension),c_new(m_c_dimension);
        compute_constraints(c0,x0);
        decision_vector x_new = x0;
        iGfun.resize(0);jGvar.resize(0); lenG=0;

        for (decision_vector::size_type j=0;j<Dc;++j)
        {
                //we perturb the component of x0 only if ub>lb, if ub=lb the variable is assumed
                //to be 'just' a parameter ... in some problem implementations this is rather
                //useful, but it also requires that the algorithm treat those variables accordingly (i.e.
                //it does not allow them to go outside the box bounds)
                if (m_ub[j] == m_lb[j]) continue;
                x_new[j] = boost::uniform_real<double>(m_lb[j],m_ub[j])(drng);
                objfun(f_new,x_new);
                compute_constraints(c_new,x_new);
                for (fitness_vector::size_type i=0;i<m_f_dimension;++i)
                {
                        if (f_new[i]!=f0[i]) {iGfun.push_back(i); jGvar.push_back(j); lenG++;}
                }
                for (constraint_vector::size_type i=0;i<m_c_dimension;++i)
                {
                        if (c_new[i]!=c0[i]) {iGfun.push_back(i+m_f_dimension); jGvar.push_back(j); lenG++;}
                }
                x_new[j] = x0[j];
        }
}


void base::set_best_x(const std::vector<decision_vector>& best_x)
{
        if(best_x.size() != 0){
                size_type n_opts = best_x.size();
                m_best_x.resize(n_opts);
                m_best_f.resize(n_opts);
                m_best_c.resize(n_opts);
                for (size_type i=0; i<n_opts; i++){
                        if(best_x.at(i).size() != get_dimension())
                                pagmo_throw(value_error,"invalid size(s) for best known decision vector(s)");
                        else{
                                //save in the class data member the value of the decision variable of solution i
                                m_best_x.at(i) = best_x.at(i);
                                //save in the class data member the corresponding value of objective(s)
                                m_best_f.at(i) = objfun(best_x.at(i));
                                //save in the class data member the corresponding value of constraint(s)
                                if(m_c_dimension>0)
                                        m_best_c.at(i) = compute_constraints(best_x.at(i));
                        }
                }
        }
}


const std::vector<constraint_vector>& base::get_best_c(void) const
{
        return m_best_c;
}


const std::vector<decision_vector>& base::get_best_x(void) const
{
        return m_best_x;
}


const std::vector<fitness_vector>& base::get_best_f(void) const
{
        return m_best_f;
}


void base::pre_evolution(population &pop) const
{
        (void)pop;
}


void base::post_evolution(population &pop) const
{
        (void)pop;
}


void base::reset_caches() const
{
        m_decision_vector_cache_f = decision_vector_cache_type(boost::numeric_cast<decision_vector_cache_type::size_type>(cache_capacity));
        m_fitness_vector_cache = fitness_vector_cache_type(boost::numeric_cast<fitness_vector_cache_type::size_type>(cache_capacity));
        m_decision_vector_cache_c = decision_vector_cache_type(boost::numeric_cast<decision_vector_cache_type::size_type>(cache_capacity));
        m_constraint_vector_cache = constraint_vector_cache_type(boost::numeric_cast<constraint_vector_cache_type::size_type>(cache_capacity));
}

}} //namespaces