cgal/QP_solver/include/CGAL/QP_solver/QP__partial_base.h

// Copyright (c) 1997-2007  ETH Zurich (Switzerland).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
// Author(s)     : Sven Schoenherr
//                 Bernd Gaertner <gaertner@inf.ethz.ch>
//                 Franz Wessendorp
//                 Kaspar Fischer

#ifndef CGAL_QP__PARTIAL_BASE_H
#define CGAL_QP__PARTIAL_BASE_H

#include <CGAL/license/QP_solver.h>


// includes
#include <CGAL/QP_solver/QP_pricing_strategy.h>
#include <CGAL/Random.h>
#include <CGAL/algorithm.h>
#include <vector>
#include <cmath>

namespace CGAL {

// ==================
// class declarations
// ==================
template < typename Q, typename ET, typename Tags >
class QP__partial_base;

template < class Solver >
struct transition_sync_functor {

    // Note that we rely here on working_vars being the number of
    // variables without artificials and the solvers in_B variable
    // being up to date. Furthermore the operator() below relies on short
    // circuit evaluation

    transition_sync_functor( const Solver& s, int w) : amb_solver(s),
        working_vars(w) { }

    bool operator() (int i) {
        return (i < working_vars) && !amb_solver.is_basic(i);
    }

private:
    const Solver& amb_solver;
    int working_vars;
};


// ===============
// class interface
// ===============
template < typename Q, typename ET, typename Tags >
class QP__partial_base : virtual public QP_pricing_strategy<Q,ET,Tags> {

    // self
    typedef  QP__partial_base<Q,ET,Tags>     Self;
    typedef  QP_pricing_strategy<Q,ET,Tags>  Base;

    typedef  typename Base::QP_solver   QP_solver;
  protected:

    // types
    typedef  typename QP_solver::Indices               Indices;
    typedef  typename QP_solver::Index_iterator        Index_iterator;
    typedef  typename QP_solver::Index_const_iterator  Index_const_iterator;

    // construction
    QP__partial_base( bool  randomize, Random&  random);

    // initialization
    virtual  void  init( );

    // access
    Index_const_iterator    active_set_begin( ) const { return N.begin();   }
    Index_const_iterator    active_set_end  ( ) const { return N.begin()+s; }

    Index_const_iterator  inactive_set_begin( ) const { return N.begin()+s; }
    Index_const_iterator  inactive_set_end  ( ) const { return N.end  ();   }

    // operations
    void  entering_basis( Index_const_iterator  it);
    void  activating    ( Index_const_iterator& it);

    virtual  void  leaving_basis( int i);
    virtual  void  transition( );

  private:

    // data members
    Indices                  N;         // non-basis;
    int                      s;         // size of active set

    bool                     permute;   // flag: permute initial non-basis
    Random&                  rand_src;  // random source

    //basic_functor<QP_solver> is_non_basic;
};

// ----------------------------------------------------------------------------

// ===============================
// class implementation (template)
// ===============================

// construction
template < typename Q, typename ET, typename Tags >  inline
QP__partial_base<Q,ET,Tags>::
QP__partial_base( bool  randomize, Random&  random)
    : permute( randomize), rand_src( random)
{ }

// initialization
template < typename Q, typename ET, typename Tags >
void
QP__partial_base<Q,ET,Tags>::
init( )
{
    // initialize indices of non-basic variables
    int  w = this->solver().number_of_working_variables();
    int  b = this->solver().number_of_basic_variables();

    if ( ! N.empty()) N.clear();
    N.reserve( w-b);                                        // use 'w' ???
    for ( int i = 0; i < w; ++i) {
        if ( ! this->solver().is_basic( i)) N.push_back( i);
    }
    if ( permute) CGAL::cpp98::random_shuffle( N.begin(), N.end(), rand_src);

    // initialize size of active set
    int  n = this->solver().number_of_variables();
    int  m = this->solver().number_of_constraints();
    // we also want to cover the high constraints/variable ratio
    if (n < m) (std::swap)(n,m);

    s = (std::min)( static_cast< unsigned int>( m*std::sqrt( n/2.0)),
                    static_cast< unsigned int>(N.size()));

    //is_non_basic.init(this->solver());
}

// operations
template < typename Q, typename ET, typename Tags >  inline
void
QP__partial_base<Q,ET,Tags>::
entering_basis( Index_const_iterator it)
{
    CGAL_qpe_precondition( it >= active_set_begin() && it < active_set_end());

    // remove from active set
    --s;
    const_cast< typename Indices::value_type&>( *it) = N[ s];
    N[ s] = N.back();
    N.pop_back();
}

template < typename Q, typename ET, typename Tags >  inline
void
QP__partial_base<Q,ET,Tags>::
activating( Index_const_iterator& it)
{
    CGAL_qpe_precondition(
        it >= inactive_set_begin() && it < inactive_set_end());

    // 'append' to active set
    std::swap( const_cast< typename Indices::value_type&>( *it), N[ s]);
    it = N.begin()+s;
    ++s;
}

template < typename Q, typename ET, typename Tags >
void
QP__partial_base<Q,ET,Tags>::
leaving_basis( int i)
{
    // all non-basic variables active?
    if ( s == static_cast< int>( N.size())) {

        // append at the end of all non-basic variables
        N.push_back( i);

    } else {

        // insert at the end of the current active subset
        N.push_back( N[ s]);
        N[ s] = i;
    }
    ++s;
}


template < typename Q, typename ET, typename Tags >
void
QP__partial_base<Q,ET,Tags>::
transition( )
{
    // Remove from N nonbasic slack and original variables that have become
    // basic during the expelling of artificial variables out of the basis
    // (between phaseI and phaseII), since these formerly nonbasic variables
    // have not entered the basis through pricing the set N has not accordingly
    // been updated.
    // Remove from N the artificial variables as well.
    // Note that we rely on the number of working variables including only
    // original and slack variables, the solvers in_B variable must be
    // up to date.
    // Furthermore we rely on std::partition not destroying the randomness
    // of the order of the nonbasic variables in N.

    int  w = this->solver().number_of_working_variables();
    transition_sync_functor<QP_solver> is_non_basic(this->solver(), w);
    N.erase( std::partition( N.begin(), N.end(),
                             is_non_basic), N.end());

    // initialize size of active set
    int  n = this->solver().number_of_variables();
    int  m = this->solver().number_of_constraints();

    s = (std::min)( static_cast< unsigned int>( m*std::sqrt( n/2.0)),
                    static_cast< unsigned int>(N.size()));
}

} //namespace CGAL

#endif // CGAL_QP__PARTIAL_BASE_H

// ===== EOF ==================================================================