minmonte.h

Go to the documentation of this file.

#ifndef MINMONTE_H
#define MINMONTE_H

#include <algorithm>
#include <cassert>
#include <iostream>

using namespace std;


// Toggle debug
#define DEBUG_MINMONTE


// Brief: Helper class: associate a value with an index.
template< typename T >
class valueindex
{
public:

  unsigned int index;
  T val;

  bool const operator < (valueindex<T> const & x)
    { return val < x.val; }

  ostream& print(ostream& os) const;

};

template< typename T >
ostream& operator << (ostream& os, valueindex<T> const & x) { return x.print(os); }


// Brief: Mimimization algorithm by an adaptive monte carlo method.
//   Assumes the minimum is in the bounding volume. 

template
< 
  typename F,   // Function functional object ()(T&,T[]).
  typename R,   // Random number functional object T r().         
  typename T,   // Number type.
  unsigned int Arraysz, // Maximum number of sampled points.
  unsigned int Dim  // Dimension (number of independent variables).
>
class minmonte
{
protected:
  

  F func;  // Function to be minimized.
  R r;     // Random number in [0,1].

  T a[Dim];  // Lower bounding box.
  T b[Dim];  // Upper bounding box.
  T lenab[Dim];  // b-a

  unsigned int s;  // The current number of samples.
  T v[Arraysz][Dim];        // Points.
  valueindex<T> pti[Arraysz];  // Value of functon applied to point.

  void calcboundingbox();

  void printv(unsigned int k);


public:

  // Constructors.
  minmonte() {}
  minmonte(F func_, R r_) : func(func_), r(r_) {}

  // Client must configure the bounding box/volume.
  template< typename V >
  void setboundingbox( V const & a_, V const & b_);
  template< typename V >
  void getboundingbox( V & a_, V & b_ ) const;

  void setfunction( F const & func_);

  // Access the top s results from the random sample. 
  //   k=0..s-1 sorted, k=0 is best.
  void getresult( T& val, T* p, unsigned int const k ) const;

  // Evaluate one iteration of the adaptive monte carlo method.
  // ns is the number of samples.
  // s is the best s cases are used to build a bounding box.
  void operator () (unsigned int s_, unsigned int ns = Arraysz );

  //
  void printab() const;

};


//
// Implementation
// ---------------------------------------------------------------

template< typename F, typename R, typename T, unsigned int Arraysz, unsigned int Dim >
void minmonte<F,R,T,Arraysz,Dim>::getresult
(
  T& val, 
  T* p,
  unsigned int const k 
) const
{
#ifdef DEBUG_MINMONTE
assert(k<s);
#endif

  val = pti[k].val;
  unsigned int index = pti[k].index;
  for (unsigned int i=0; i<Dim; ++i, ++p)    
    *p = v[index][i]; 
}
  


template< typename F, typename R, typename T, unsigned int Arraysz, unsigned int Dim >
template< typename V >
void minmonte<F,R,T,Arraysz,Dim>::getboundingbox(V & a_, V & b_) const
{
  for (unsigned int i=0; i<Dim; ++i)
  {
    a_[i] = a[i];
    b_[i] = b[i];
  }
}


template< typename F, typename R, typename T, unsigned int Arraysz, unsigned int Dim >
template< typename V >
void minmonte<F,R,T,Arraysz,Dim>::setboundingbox(V const & a_, V const & b_)
{
  for (unsigned int i=0; i<Dim; ++i)
  {
    a[i] = a_[i];
    b[i] = b_[i];
    lenab[i] = b[i]-a[i];
#ifdef DEBUG_MINMONTE
assert(a[i]<b[i]);
#endif
  }

}

template< typename F, typename R, typename T, unsigned int Arraysz, unsigned int Dim >
void minmonte<F,R,T,Arraysz,Dim>::setfunction(F const & func_)
{
  func = func_;
} 


template< typename F, typename R, typename T, unsigned int Arraysz, unsigned int Dim >
void minmonte<F,R,T,Arraysz,Dim>::operator () (unsigned int s_, unsigned int ns )
{
  
#ifdef DEBUG_MINMONTE

if (false)
  cout << "Entering minmonte ()" << endl;

assert( ns <= Arraysz );
assert(s_>0);
assert(s_<=ns);

#endif

  s = s_;

  //<TODO> I had trouble with pointer arithmetic.
  // Possibly the use of inline function? In the mean time pass 
  //   a single array.  For some reason passing a double array as a single failed.
  T temp[Dim];

  for (unsigned int i=0; i<ns; ++i)
  {
    for (unsigned int k=0; k<Dim; ++k)
      temp[k] = v[i][k] = a[k] + lenab[k]*r();

    func(pti[i].val,temp);
    //func(pti[i].val,(T*)(v+i*Dim));
    pti[i].index = i;

#ifdef DEBUG_MINMONTE
if (false)
{
cout << pti[i];    
printv(i);
cout << endl;
}
#endif

  };

  partial_sort(pti,pti+s,pti+ns);

    
#ifdef DEBUG_MINMONTE
if (false)
{
  cout << "After partial sort." << endl;
  for (unsigned int i=0; i<s; ++i)
  {
    cout << "pti[" << i << "]=" << pti[i] << " ";
    printv(pti[i].index);
    cout << endl;
  }
}
#endif

  calcboundingbox();

}


template< typename F, typename R, typename T, unsigned int Arraysz, unsigned int Dim >
void minmonte<F,R,T,Arraysz,Dim>::calcboundingbox()
{
  // Assume a and b are valid.
  for (unsigned int k=0; k<Dim; ++k)
    a[k] = b[k] = v[pti[0].index][k];

#ifdef DEBUG_MINMONTE
if(false)
{
  cout << "***" << endl;
  printv( pti[0].index );
  cout << endl;
  printab();
}
#endif

  for (unsigned int k=0; k<Dim; ++k)
  {
    for (unsigned int i=1; i<s; ++i)
    {
      a[k] = min(a[k],v[pti[i].index][k]);
      b[k] = max(b[k],v[pti[i].index][k]);
    }
  }

  for (unsigned int k=0; k<Dim; ++k)
    lenab[k] = b[k]-a[k];
}


template< typename F, typename R, typename T, unsigned int Arraysz, unsigned int Dim >
void minmonte<F,R,T,Arraysz,Dim>::printv(unsigned int k)
{
  cout << "v[" << k << "]=[";
  for (unsigned int i=0; i<Dim-1; ++i)
        cout << v[k][i] << ",";
  cout << v[k][Dim-1] << "]";
}

template< typename F, typename R, typename T, unsigned int Arraysz, unsigned int Dim >
void minmonte<F,R,T,Arraysz,Dim>::printab() const
{
  cout << "a=[";
  for (unsigned int i=0; i<Dim-1; ++i)
    cout << a[i] << ",";
  cout << a[Dim-1] << "]" << "  ";
  cout << "b=[";
  for (unsigned int i=0; i<Dim-1; ++i)
    cout << b[i] << ",";
  cout << b[Dim-1] << "]" << "  ";
  cout << "lenab=[";
  for (unsigned int i=0; i<Dim-1; ++i)
    cout << lenab[i] << ",";
  cout << lenab[Dim-1] << "]" << endl;
}

template< typename T >
ostream& valueindex<T>::print(ostream& os) const
{
  os << "(" << index << "," << val << ")";
  return os;
}



#endif