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

4 files changed, 364 insertions, 0 deletions

diff --git a/methods/linmin/bas_util.hpp b/methods/linmin/bas_util.hpp
new file mode 100644
index 0000000..ccdce9e
--- /dev/null
+++ b/methods/linmin/bas_util.hpp
@@ -0,0 +1,65 @@
+#ifndef BAS_UTIL
+#define BAS_UTIL
+#define OPT_HEADER
+#include <core/opt_traits.hpp>
+#include <algorithm>
+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/linmin/brent.hpp b/methods/linmin/brent.hpp
new file mode 100644
index 0000000..1a0ee39
--- /dev/null
+++ b/methods/linmin/brent.hpp
@@ -0,0 +1,111 @@
+#ifndef BRENT_HPP
+#define BRENT_HPP
+#define OPT_HEADER
+#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/linmin/linmin.hpp b/methods/linmin/linmin.hpp
new file mode 100644
index 0000000..ddf2a8e
--- /dev/null
+++ b/methods/linmin/linmin.hpp
@@ -0,0 +1,106 @@
+#ifndef LINMIN_HPP
+#define LINMIN_HPP
+#define OPT_HEADER
+#include "mnbrak.hpp"
+#include "brent.hpp"
+#include <cmath>
+#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&)
+      :func_obj<rT,rT>(),p1(),xi1(),pfoo(0)
+    {}
+    
+  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=std::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/linmin/mnbrak.hpp b/methods/linmin/mnbrak.hpp
new file mode 100644
index 0000000..9871473
--- /dev/null
+++ b/methods/linmin/mnbrak.hpp
@@ -0,0 +1,82 @@
+#ifndef MNBRAK_HPP
+#define MNBRAK_HPP 
+#define OPT_HEADER
+//#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