---

total_energy_function.cpp

//$Header: /home/ben/Mapper/c++/RCS/total_energy_function.cpp,v 6.4 2002/06/25 16:13:50 ben Exp $
// Copyright Benedict Adamson 2002.
// This file is part of Mapper.

// Mapper 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.

// Mapper 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 Mapper; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

#include <climits>

#include "total_energy_function.h"



const unsigned total_energy_function::fixed_point = UINT_MAX;



total_energy_function::total_energy_function(
   const double t
   )
   :
   tol_( t ),
   points_( 0 ),
   fixed_points_( 0 ),
   constraints_( 0 ),
   constraint_to_var(),
   point_to_var()
{
   {//preconditions:
      assert( 0 < t && t < 1 );
   };
   {//postconditions
      assert( this->tol() == t );
      assert( this->points() == 0 );
      assert( this->fixed_points() == 0 );
      assert( this->constraints() == 0 );
   };
}



total_energy_function::~total_energy_function(void)
{
   //Do nothing
}



void total_energy_function::set_tol(
   const double t
   )
{
   {//preconditions:
      assert( 0 < t && t < 1 );
   };
   this->tol_ = t;
   {//preconditions:
      assert( this->tol() == t );
   };
}



void assoc(
   total_energy_function &f,
   const vector< vector2 * > &po,
   const set< vector2 * > &fp,
   const vector< constraint * > &co
   )
{
   {//preconditions
      assert( f.points() == 0 );
      assert( f.fixed_points() == 0 );
      assert( f.constraints() == 0 );
      assert( fp.size() <= po.size() );
   };
   f.points_ = &po;
   f.fixed_points_ = &fp;
   f.constraints_ = &co;
   {//compute point_to_var
      unsigned i = 0;
      for( vector< vector2 * >::const_iterator j = po.begin();
           j != po.end();
           j++ )
      {
         vector2 *pj = *j;
         assert( pj );
         if( fp.count( pj ) ) //pj is a fixed point
            f.point_to_var[ pj ] = total_energy_function::fixed_point;
         else{
            f.point_to_var[ pj ] = i;
            i += 2; //allocate variables for x and y coordinates
         };
      };
   };
   {//compute constraint_to_var
      const unsigned nc = co.size();
      f.constraint_to_var = vector< vector<unsigned> >( nc );
      for( unsigned i = 0; i < nc; i++ ){
         const vector< vector2 * > &cp( co[i]->points() );
         const unsigned np = cp.size();
         assert( 0 < np ); //otherwise, constrains nothing
         f.constraint_to_var[i] = vector<unsigned>( np );
         for( unsigned j = 0; j < np; j++ ){
            vector2 *p = cp[j];
            assert( p );
            assert( f.point_to_var.count( p ) );
            f.constraint_to_var[i][j] = f.point_to_var[ p ];
         };
      };
   };
   {//postconditions
      assert( f.points() == &po );
      assert( f.fixed_points() == &fp );
      assert( f.constraints() == &co );
   };
}



double total_energy_function::operator() (
   const vector<double> &x
   ) const
{
   {//preconditions:
      assert( this->points() );
      assert( this->fixed_points() );
      assert( this->constraints() );
      assert( x.size() == this->n_vars() );
   };
   double e;
   vector<double> dedx( this->n_vars() );
   this->evaluate( x, e, dedx );
   //discard dedx
   return e;
}




void total_energy_function::evaluate(
   const vector<double> &x,
   double &e, //output
   vector<double> &gradient //output
   ) const
{
   {//preconditions:
      assert( this->points() );
      assert( this->fixed_points() );
      assert( this->constraints() );
      assert( x.size() == this->n_vars() );
   };
   //f == (*this)( x )
   {//initialise running totals
      e = 0.0;
      gradient = vector<double>( this->n_vars() );
   };
   const unsigned n_constraints = this->constraints()->size();
   const double tol_c = this->tol()/double(1U + n_constraints);
   for( unsigned i = 0; i < n_constraints; i++ ){
      const constraint *c =  (*(this->constraints()))[i];
      assert( c );
      const vector< unsigned > &ctv( this->constraint_to_var[i] );
      const unsigned np = c->points().size();
      assert( 0 < np );
      vector< vector2 > p( np ), dedp( np );
      for( unsigned j = 0; j < np; j++ ){
         const unsigned var = ctv[j];
         if( var == fixed_point ){
            const vector2 *pj = c->points()[j];
            assert( pj );
            p[j] = *pj;
         }else{
            assert( var+1 < this->n_vars() );
            p[j] = vector2( x[ var ], x[ var+1 ] );
         };
      };
      double ep;
      c->energy( tol_c, p, ep, dedp );
      {//increase the running totals
         e += ep;
         for( unsigned j = 0; j < np; j++ ){
            const unsigned var = ctv[j];
            if( var != fixed_point ){
               assert( var+1 < this->n_vars() );
               gradient[ var ] += dedp[j].x;
               gradient[ var+1 ] += dedp[j].y;
            };//else no effect on gradient
         };
      };
   };
}



vector<double> total_energy_function::current_var(void) const
{
   {//preconditions
      assert( this->points() );
      assert( this->fixed_points() );
   };
   vector<double> result( this->n_vars() );
   map< vector2 *, unsigned > &m(
      const_cast< map< vector2 *, unsigned > & >(
         this->point_to_var
         )); //circumvent const
   for( vector< vector2 *const >::iterator i = points()->begin();
        i != points()->end();
        i++ ){
         vector2 *p = *i;
         assert( p );
         assert( m.count( p ) );
         const unsigned j = m[ p ];
         if( j == fixed_point )
            ;//no variables
         else{ //j and j+1 are the corresponding variables
            assert( j+1 < this->n_vars() );
            result[j] = p->x;
            result[j+1] = p->y;
         };
   };
   return result;
}



void total_energy_function::move_points(
   const vector<double> & x
   ) const
{
   {//preconditions:
      assert( x.size() == this->n_vars() );
   };
   map< vector2 *, unsigned > &m(
      const_cast< map< vector2 *, unsigned > & >(
         this->point_to_var
         )); //circumvent const
   for( vector< vector2 *const >::iterator i = points()->begin();
        i != points()->end();
        i++ ){
         vector2 *p = *i;
         assert( p );
         assert( m.count( p ) );
         const unsigned j = m[ p ];
         if( j == fixed_point )
            ;//no variables
         else{ //j and j+1 are the corresponding variables
            assert( j+1 < this->n_vars() );
            p->x = x[j];
            p->y = x[j+1];
         };
   };
   {//postconditions:
      assert( this->current_var() == x );
   };
}

---