git-svn-id: file:///home/svn/opt_utilities@1 ed2142bd-67ad-457f-ba7c-d818d4011675

6 files changed, 806 insertions, 0 deletions

diff --git a/methods/gsl_simplex/gsl_simplex.hpp b/methods/gsl_simplex/gsl_simplex.hpp
new file mode 100644
index 0000000..442c9f0
--- /dev/null
+++ b/methods/gsl_simplex/gsl_simplex.hpp
@@ -0,0 +1,199 @@
+#ifndef GSL_SIMPLEX_METHOD
+#define GSL_SIMPLEX_METHOD
+#include <core/optimizer.hpp>
+//#include <blitz/array.h>
+#include <vector>
+#include <limits>
+#include <cassert>
+#include <cmath>
+#include <algorithm>
+#include <gsl_multimin.h>
+/*
+ *
+*/
+#include <iostream>
+
+
+namespace opt_utilities
+{
+  template <typename rT,typename pT>
+  double gsl_func_adapter(const gsl_vector* v,void* params)
+  {
+    pT temp;
+    temp.resize(v->size);
+    for(size_t i=0;i<get_size(temp);++i)
+      {
+	set_element(temp,i,gsl_vector_get(v,i));
+      }
+    return ((func_obj<rT,pT>*)params)->eval(temp);
+  }
+
+
+  template <typename rT,typename pT>
+  class gsl_simplex
+    :public opt_method<rT,pT>
+  {
+  public:
+    typedef pT array1d_type;
+    typedef rT T;
+  private:
+    func_obj<rT,pT>* p_fo;
+    optimizer<rT,pT>* p_optimizer;
+    
+    //typedef blitz::Array<rT,2> array2d_type;
+    
+    
+  private:
+    array1d_type start_point;
+    array1d_type end_point;
+    
+  private:
+    rT threshold;
+  private:
+    rT func(const pT& x)
+    {
+      assert(p_fo!=0);
+      return p_fo->eval(x);
+    }
+
+   
+  public:
+    gsl_simplex()
+      :threshold(1e-4)
+    {}
+
+    virtual ~gsl_simplex()
+    {
+    };
+
+    gsl_simplex(const gsl_simplex<rT,pT>& rhs)
+      :p_fo(rhs.p_fo),p_optimizer(rhs.p_optimizer),
+       start_point(rhs.start_point),
+       end_point(rhs.end_point),
+       threshold(rhs.threshold)
+    {
+    }
+
+    gsl_simplex<rT,pT>& operator=(const gsl_simplex<rT,pT>& rhs)
+    {
+      threshold=rhs.threshold;
+      p_fo=rhs.p_fo;
+      p_optimizer=rhs.p_optimizer;
+      opt_eq(start_point,rhs.start_point);
+      opt_eq(end_point,rhs.end_point);
+    }
+    
+    opt_method<rT,pT>* do_clone()const
+    {
+      return new gsl_simplex<rT,pT>(*this);
+    }
+    
+    void do_set_start_point(const array1d_type& p)
+    {
+      start_point.resize(get_size(p));
+      opt_eq(start_point,p);
+      
+    }
+
+    void do_set_precision(rT t)
+    {
+      threshold=t;
+    }
+
+    void do_set_optimizer(optimizer<rT,pT>& o)
+    {
+      p_optimizer=&o;
+      p_fo=p_optimizer->ptr_func_obj();
+    }
+    
+    
+    
+    pT do_optimize()
+    {
+      const gsl_multimin_fminimizer_type *T = 
+	gsl_multimin_fminimizer_nmsimplex;
+      gsl_multimin_fminimizer *s = NULL;
+      gsl_vector *ss, *x;
+      gsl_multimin_function minex_func;
+      
+      size_t iter = 0;
+      int status;
+      double size;
+      
+      /* Starting point */
+      x = gsl_vector_alloc (get_size(start_point));
+      //      gsl_vector_set (x, 0, 5.0);
+      //gsl_vector_set (x, 1, 7.0);
+      for(size_t i=0;i!=get_size(start_point);++i)
+	{
+	  gsl_vector_set(x,i,get_element(start_point,i));
+	}
+
+
+      /* Set initial step sizes to 1 */
+      ss = gsl_vector_alloc (get_size(start_point));
+      gsl_vector_set_all (ss, 1.0);
+
+
+      //foo f;
+      /* Initialize method and iterate */
+      minex_func.n = get_size(start_point);
+      minex_func.f = &gsl_func_adapter<double,std::vector<double> >;
+      minex_func.params = (void *)p_fo;
+      
+      s = gsl_multimin_fminimizer_alloc (T, get_size(start_point));
+      gsl_multimin_fminimizer_set (s, &minex_func, x, ss);
+      
+      do
+	{
+	  iter++;
+	  status = gsl_multimin_fminimizer_iterate(s);
+	  
+	  if (status) 
+	    break;
+	  //std::cerr<<"threshold="<<threshold<<std::endl;
+	  size = gsl_multimin_fminimizer_size (s);
+	  status = gsl_multimin_test_size (size, threshold);
+	  
+	  if (status == GSL_SUCCESS)
+	    {
+	      //printf ("converged to minimum at\n");
+	    }
+	  
+	  //printf ("%5d %10.3e %10.3ef f() = %7.3f size = %.3f\n", 
+	  //iter,
+	  //gsl_vector_get (s->x, 0), 
+	  //gsl_vector_get (s->x, 1), 
+	  //  s->fval, size);
+	}
+      while (status == GSL_CONTINUE && iter < 100);
+      
+      /*
+	foo f;
+	gsl_vector_set (x, 0, 0.0);
+	gsl_vector_set (x, 1, 0.0);
+	cout<<"fdsa ";
+	cout<<gsl_func_adapter<double,vector<double> >(x,(void*)&f)<<endl;;
+	
+      */
+      
+      end_point.resize(get_size(start_point));
+      for(size_t i=0;i<get_size(start_point);++i)
+	{
+	  set_element(end_point,i,gsl_vector_get(s->x,i));
+	}
+
+      gsl_vector_free(x);
+      gsl_vector_free(ss);
+      gsl_multimin_fminimizer_free (s);
+      
+      
+      return end_point;
+    } 
+  };
+  
+};
+
+
+#endif
+//EOF
diff --git a/methods/powell/bas_util.hpp b/methods/powell/bas_util.hpp
new file mode 100644
index 0000000..743dc3c
--- /dev/null
+++ b/methods/powell/bas_util.hpp
@@ -0,0 +1,63 @@
+#ifndef BAS_UTIL
+#define BAS_UTIL
+#include <core/opt_traits.hpp>
+namespace opt_utilities
+{
+  template <typename T>
+  T tabs(T x)
+  {
+    return T(x)<T(0)?T(-x):T(x);
+  }
+
+  template <typename T>
+  T sqr(T x)
+  {
+    return x*x;
+  }
+  
+  
+  template <typename T>
+  void shft3(T&a,T& b,T& c,T d)
+  {
+    opt_eq(a,b);
+    opt_eq(b,c);
+    opt_eq(c,d);
+  }
+
+  template <typename T>
+  void shft(T& a,T& b,T& c,T d)
+  {
+    opt_eq(a,b);
+    opt_eq(b,c);
+    opt_eq(c,d);
+  }
+  template <typename T>  
+  void swap(T& ax,T& bx)
+  {
+    //  swap(ax,bx);
+    T temp;
+    opt_eq(temp,ax);
+    opt_eq(ax,bx);
+    opt_eq(bx=temp);
+  }
+  
+  template <typename T>
+  T sign(const T& a,const T& b)
+  {
+    return b>=0?T(a>=0?T(a):T(-a)):T(a>=0?T(-a):T(a));
+  }
+  
+  template <typename T>
+  T max(T a,T b)
+  {
+    return b>a?T(b):T(a);
+  }
+
+  template <typename T>
+  void mov3(T& a,T& b,T& c, T& d,T& e,T& f)
+  {
+    opt_eq(a,d);opt_eq(b,e);opt_eq(c,f);
+  }
+}
+
+#endif
diff --git a/methods/powell/brent.hpp b/methods/powell/brent.hpp
new file mode 100644
index 0000000..cb3b962
--- /dev/null
+++ b/methods/powell/brent.hpp
@@ -0,0 +1,110 @@
+#ifndef BRENT_HPP
+#define BRENT_HPP
+#include <iostream>
+#include "bas_util.hpp"
+//#include "optimizer.hpp"
+namespace opt_utilities
+{
+  template<typename T>
+  T brent(T ax,T bx,T cx,func_obj<T,T>& f,T tol,T& xmin)
+  {
+    const int ITMAX=100;
+    const T CGOLD=0.3819660;
+    const T ZEPS=std::numeric_limits<T>::epsilon()*1.e-3;
+    
+    int iter;
+    T a=0,b=0,d(0),etemp=0,fu=0,fv=0,fw=0,fx=0,p=0,q=0
+      ,r=0,tol1=0,tol2=0,u=0,v=0,w=0,x=0,xm=0;
+    T e=0.;
+    a=(ax<cx?ax:cx);
+    b=(ax>cx?ax:cx);
+    x=w=v=bx;
+    fw=fv=fx=f.eval(x);
+    for(iter=0;iter<ITMAX;++iter)
+      {
+	xm=.5*(a+b);
+	tol2=2.*(tol1=tol*tabs(x)+ZEPS);
+	if(tabs(T(x-xm))<=(tol2-.5*(b-a)))
+	  {
+	    xmin=x;
+	    return fx;
+	  }
+	if(tabs(e)>tol1)
+	  {
+	    r=(x-w)*(fx-fv);
+	    q=(x-v)*(fx-fw);
+	    p=(x-v)*q-(x-w)*r;
+	    q=2.*(q-r);
+	    if(q>0.)
+	      {
+		p=-p;
+	      }
+	    q=tabs(q);
+	    etemp=e;
+	    e=d;
+	    if(tabs(p)>=tabs(T(T(.5)*p*etemp))||p<=q*(a-x)||p>=q*(b-x))
+	      {
+		d=CGOLD*(e=(x>=xm?a-x:b-x));
+	      }
+	    else
+	      {
+		d=p/q;
+		u=x+d;
+		if(u-a<tol2||b-u<tol2)
+		  {
+		    d=sign(tol1,T(xm-x));
+		  }
+	      }
+	    
+	  }
+	else
+	  {
+	    d=CGOLD*(e=(x>=xm?a-x:b-x));
+	  }
+	u=(tabs(d)>=tol1?x+d:x+sign(tol1,d));
+	fu=f.eval(u);
+	if(fu<=fx)
+	  {
+	    if(u>=x)
+	      {
+		a=x;
+	      }
+	    else
+	      {
+		b=x;
+	      }
+	    shft3(v,w,x,u);
+	    shft3(fv,fw,fx,fu);
+	  }
+	else
+	  {
+	    if(u<x)
+	      {
+		a=u;
+	      }
+	    else
+	      {
+		b=u;
+	      }
+	    if(fu<=fw||w==x)
+	      {
+		v=w;
+		w=u;
+		fv=fw;
+		fw=fu;
+	      }
+	    else if(fu<=fv||v==x||v==w)
+	      {
+		v=u;
+		fv=fu;
+	      }
+	  }
+      }
+    std::cerr<<"Too many iterations in brent"<<std::endl;
+    xmin=x;
+    return fx;
+    
+  }
+}
+
+#endif
diff --git a/methods/powell/linmin.hpp b/methods/powell/linmin.hpp
new file mode 100644
index 0000000..5ab5af0
--- /dev/null
+++ b/methods/powell/linmin.hpp
@@ -0,0 +1,102 @@
+#ifndef LINMIN_HPP
+#define LINMIN_HPP
+#include "mnbrak.hpp"
+#include "brent.hpp"
+#include <core/opt_traits.hpp>
+
+namespace opt_utilities
+{
+  template <typename rT,typename pT>
+  class func_adaptor
+    :public func_obj<rT,rT>
+  {
+  private:
+    const pT p1,xi1;
+    const func_obj<rT,pT>* pfoo;
+    func_adaptor(){}
+    func_adaptor(const func_adaptor&){}
+    
+  public:
+    /*
+    void set_origin(pT& p2)
+    {
+      p1=p2;
+    }
+    
+    void set_direction(pT& xi2)
+    {
+      xi1=xi2;
+    }
+    
+    void set_func_obj(func_obj<rT,pT>& foo)
+    {
+      pfoo=&foo;
+      }*/
+  public:
+    func_adaptor(const pT& _p,const pT& _xi,const func_obj<rT,pT>& pf)
+      :p1(_p),xi1(_xi),pfoo(&pf)
+    {}
+    
+  private:
+    func_obj<rT,rT>* do_clone()const
+    {
+      return new func_adaptor(*this);
+    }
+      
+    rT do_eval(const rT& x)
+    {
+      //assert(p1.size()==xi1.size());
+      
+      pT xt;
+      opt_eq(xt,p1);
+      for(size_t i=0;i<get_size(xt);++i)
+	{
+	  //get_element(xt,i)+=x*get_element((pT)xi1,i);
+	  set_element(xt,i,
+		      get_element(xt,i)+x*get_element((pT)xi1,i));
+	  //get_element((pT)xi1,i);
+	}
+      return const_cast<func_obj<rT,pT>&>(*pfoo).eval(xt);
+      //return x;
+    }
+  };
+  
+
+  template<typename rT,typename pT>
+  void linmin(pT& p,pT& xi,rT& fret,func_obj<rT,pT>& func)
+  {
+
+    //  assert(p.size()==10);
+    //assert(xi.size()==10);
+    func_adaptor<rT,pT> fadpt(p,xi,func);
+
+    int j=0;
+    const rT TOL=sqrt(std::numeric_limits<rT>::epsilon());
+    rT xx=0,xmin=0,fx=0,fb=0,fa=0,bx=0,ax=0;
+    int n=(int)get_size(p);
+
+
+    ax=0.;
+    xx=1.;
+    
+
+    mnbrak(ax,xx,bx,fa,fx,fb,fadpt);
+    //cout<<xx<<endl;
+    fret=brent(ax,xx,bx,fadpt,TOL,xmin);
+    //cout<<xmin<<endl;
+    for(j=0;j<n;++j)
+      {
+	//get_element(xi,j)*=xmin;
+	set_element(xi,j,
+		    get_element(xi,j)*xmin);
+	//get_element(p,j)+=get_element(xi,j);
+	set_element(p,j,
+		    get_element(p,j)+get_element(xi,j));
+      }
+    //  delete xicom_p;
+    //delete pcom_p;
+  }
+}
+
+
+#endif
diff --git a/methods/powell/mnbrak.hpp b/methods/powell/mnbrak.hpp
new file mode 100644
index 0000000..4887bec
--- /dev/null
+++ b/methods/powell/mnbrak.hpp
@@ -0,0 +1,81 @@
+#ifndef MNBRAK_HPP
+#define MNBRAK_HPP 
+//#include "optimizer.hpp"
+#include "bas_util.hpp"
+namespace opt_utilities
+{
+
+  
+  template <typename T>
+  void mnbrak(T& ax,T& bx,T& cx,T& fa,T& fb,T& fc,func_obj<T,T>& func)
+  {
+    const T GOLD=1.618034;
+    const T GLIMIT=100;
+    const T TINY=std::numeric_limits<T>::epsilon();
+    T ulim,u,r,q,fu;
+    fa=func.eval(ax);
+    fb=func.eval(bx);
+    
+    if(fb>fa)
+      {
+	//shft(dum,ax,bx,dum);
+	//shft(dum,fb,fa,dum);
+	std::swap(ax,bx);
+	std::swap(fa,fb);
+      }
+
+    cx=bx+GOLD*(bx-ax);
+    fc=func.eval(cx);
+    while(fb>fc)
+      {
+	r=(bx-ax)*(fb-fc);
+	q=(bx-cx)*(fb-fa);
+	u=bx-T((bx-cx)*q-(bx-ax)*r)/
+	  T(T(2.)*sign(T(max(T(tabs(T(q-r))),T(TINY))),T(q-r)));
+	ulim=bx+GLIMIT*(cx-bx);
+	if((bx-u)*(u-cx)>0.)
+	  {
+	    fu=func.eval(u);
+	    if(fu<fc)
+	      {
+		ax=bx;
+		bx=u;
+		fa=fb;
+		fb=fu;
+		return;
+	      }
+	    else if(fu>fb)
+	      {
+		cx=u;
+		fc=fu;
+		return;
+	      }
+	    u=cx+GOLD*(cx-bx);
+	    fu=func.eval(u);
+	  }
+	else if((cx-u)*(u-ulim)>0.)
+	  {
+	    fu=func.eval(u);
+	    if(fu<fc)
+	      {
+		shft3(bx,cx,u,T(cx+GOLD*(cx-bx)));
+		shft3(fb,fc,fu,func.eval(u));
+	      }
+	  }
+	else if((u-ulim)*(ulim-cx)>=0)
+	  {
+	    u=ulim;
+	    fu=func.eval(u);
+	  }
+	else
+	  {
+	    u=cx+GOLD*(cx-bx);
+	    fu=func.eval(u);
+	  }
+	shft3(ax,bx,cx,u);
+	shft3(fa,fb,fc,fu);
+      }
+  }
+}
+
+#endif
diff --git a/methods/powell/powell_method.hpp b/methods/powell/powell_method.hpp
new file mode 100644
index 0000000..bec12b7
--- /dev/null
+++ b/methods/powell/powell_method.hpp
@@ -0,0 +1,251 @@
+#ifndef POWELL_METHOD
+#define POWELL_METHOD
+#include <core/optimizer.hpp>
+//#include <blitz/array.h>
+#include <limits>
+#include <cassert>
+#include <cmath>
+#include "linmin.hpp"
+#include <algorithm>
+/*
+ *
+*/
+#include <iostream>
+
+namespace opt_utilities
+{
+  /*
+  template <typename T>
+  T tabs(T x)
+  {
+    return x<0?-x:x;
+  }
+  */
+
+  template <typename rT,typename pT>
+  class powell_method
+    :public opt_method<rT,pT>
+  {
+  public:
+    typedef pT array1d_type;
+    typedef rT T;
+  private:
+    func_obj<rT,pT>* p_fo;
+    optimizer<rT,pT>* p_optimizer;
+    
+    //typedef blitz::Array<rT,2> array2d_type;
+    
+    
+  private:
+    array1d_type start_point;
+    array1d_type end_point;
+    
+  private:
+    int ncom;
+    array1d_type pcom_p;
+    array1d_type xicom_p;
+    rT threshold;
+    T** xi;
+    T* xi_1d;
+  private:
+    rT func(const pT& x)
+    {
+      assert(p_fo!=0);
+      return p_fo->eval(x);
+    }
+
+   
+  private:
+    void clear_xi()
+    {
+      if(xi_1d!=0)
+	{
+	  delete[] xi_1d;
+	}
+      if(xi!=0)
+	{
+	  delete[] xi;
+	}
+    }
+
+    void init_xi(int n)
+    {
+      clear_xi();
+      xi_1d=new T[n*n];
+      xi=new T*[n];
+      for(int i=0;i!=n;++i)
+	{
+	  xi[i]=xi_1d+i*n;
+	}
+      for(int i=0;i!=n;++i)
+	{
+	  for(int j=0;j!=n;++j)
+	    {
+	      xi[i][j]=(j==i?1:0);
+	    }
+	}
+    }
+
+
+
+    void powell(array1d_type& p,const T ftol,
+	   int& iter,T& fret)
+    {
+      const int ITMAX=200;
+      const T TINY=std::numeric_limits<T>::epsilon();
+      int i,j,ibig;
+      T del,fp,fptt,t;
+      int n=(int)get_size(p);
+      array1d_type pt(n);
+      array1d_type ptt(n);
+      array1d_type xit(n);
+      fret=p_fo->eval(p);
+      
+      for(j=0;j<n;++j)
+	{
+	  //get_element(pt,j)=get_element(p,j);
+	  set_element(pt,j,get_element(p,j));
+	}
+      for(iter=0;;++iter)
+	{
+	  fp=fret;
+	  ibig=0;
+	  del=0.0;
+	  for(i=0;i<n;++i)
+	    {
+	      for(j=0;j<n;++j)
+		{
+		  //get_element(xit,j)=xi[j][i];
+		  set_element(xit,j,xi[j][i]);
+		}
+	      fptt=fret;
+	      linmin(p,xit,fret,(*p_fo));
+	      if((fptt-fret)>del)
+		{
+		  del=fptt-fret;
+		  ibig=i+1;
+		}
+	    }
+	  if(T(2.)*(fp-fret)<=ftol*(tabs(fp)+tabs(fret))+TINY)
+	    {
+	      return;
+	    }
+	  if(iter==ITMAX)
+	    {
+	      std::cerr<<"powell exceeding maximun iterations."<<std::endl;
+	      return;
+	    }
+	  for(j=0;j<n;++j)
+	    {
+	      //get_element(ptt,j)=T(2.)*get_element(p,j)-get_element(pt,j);
+	      set_element(ptt,j,T(2.)*get_element(p,j)-get_element(pt,j));
+	      //get_element(xit,j)=
+	      //get_element(p,j)-get_element(pt,j);
+	      set_element(xit,j,get_element(p,j)-get_element(pt,j));
+	      //get_element(pt,j)=get_element(p,j);
+	      set_element(pt,j,get_element(p,j));
+	    }
+	  fptt=func(ptt);
+	  if(fptt<fp)
+	    {
+	      t=T(2.)*(fp-T(2.)*fret+fptt)*sqr(T(fp-fret-del))-del*sqr(T(fp-fptt));
+	      if(t<T(0.))
+		{
+		  linmin(p,xit,fret,*p_fo);
+		  for(j=0;j<n;++j)
+		    {
+		      xi[j][ibig-1]=xi[j][n-1];
+		      xi[j][n-1]=get_element(xit,j);
+		      
+		    }
+		}
+	    }
+	}
+    }
+
+    
+  public:
+    
+    powell_method()
+      :threshold(1e-4),xi(0),xi_1d(0)
+    {}
+
+    virtual ~powell_method()
+    {
+      clear_xi();
+    };
+
+    powell_method(const powell_method<rT,pT>& rhs)
+      :p_fo(rhs.p_fo),p_optimizer(rhs.p_optimizer),
+       start_point(rhs.start_point),
+       end_point(rhs.end_point),
+       ncom(rhs.ncom),
+       threshold(rhs.threshold),xi(0),xi_1d(0)
+    {
+    }
+
+    powell_method<rT,pT>& operator=(const powell_method<rT,pT>& rhs)
+    {
+      threshold=rhs.threshold;
+      xi=0;
+      xi_1d=0;
+      p_fo=rhs.p_fo;
+      p_optimizer=rhs.p_optimizer;
+      start_point=rhs.start_point;
+      end_point=rhs.end_point;
+      ncom=rhs.ncom;
+      threshold=rhs.threshold;
+    }
+    
+    opt_method<rT,pT>* do_clone()const
+    {
+      return new powell_method<rT,pT>(*this);
+    }
+    
+    void do_set_start_point(const array1d_type& p)
+    {
+      resize(start_point,get_size(p));
+      opt_eq(start_point,p);
+      
+    }
+
+    void do_set_precision(rT t)
+    {
+      threshold=t;
+    }
+
+    void do_set_optimizer(optimizer<rT,pT>& o)
+    {
+      p_optimizer=&o;
+      p_fo=p_optimizer->ptr_func_obj();
+    }
+    
+    
+    
+    pT do_optimize()
+    {
+
+      init_xi((int)get_size(start_point));
+
+
+      for(int i=0;i<(int)get_size(start_point);++i)
+	{
+	  for(int j=0;j<(int)get_size(start_point);++j)
+	    {
+	      xi[i][j]=(i==j)?1:0;
+	    }
+	}
+      
+      int iter=100;
+      opt_eq(end_point,start_point);
+      rT fret;
+      powell(end_point,threshold,iter,fret);
+      return end_point;
+    } 
+  };
+
+};
+
+
+#endif
+//EOF