From 109473cacf85c7e2a2831d56531399ab4621654e Mon Sep 17 00:00:00 2001
From: astrojhgu <astrojhgu@ed2142bd-67ad-457f-ba7c-d818d4011675>
Date: Thu, 3 Jun 2010 17:55:42 +0000
Subject: git-svn-id: file:///home/svn/opt_utilities@118
 ed2142bd-67ad-457f-ba7c-d818d4011675

---
 methods/bfgs/bfgs.hpp            | 200 +++++++++++++++++++++++++++++++++++++++
 methods/linmin/bas_util.hpp      |  65 +++++++++++++
 methods/linmin/brent.hpp         | 111 ++++++++++++++++++++++
 methods/linmin/linmin.hpp        | 106 +++++++++++++++++++++
 methods/linmin/mnbrak.hpp        |  82 ++++++++++++++++
 methods/powell/bas_util.hpp      |  65 -------------
 methods/powell/brent.hpp         | 111 ----------------------
 methods/powell/linmin.hpp        | 105 --------------------
 methods/powell/mnbrak.hpp        |  82 ----------------
 methods/powell/powell_method.hpp |   2 +-
 10 files changed, 565 insertions(+), 364 deletions(-)
 create mode 100644 methods/bfgs/bfgs.hpp
 create mode 100644 methods/linmin/bas_util.hpp
 create mode 100644 methods/linmin/brent.hpp
 create mode 100644 methods/linmin/linmin.hpp
 create mode 100644 methods/linmin/mnbrak.hpp
 delete mode 100644 methods/powell/bas_util.hpp
 delete mode 100644 methods/powell/brent.hpp
 delete mode 100644 methods/powell/linmin.hpp
 delete mode 100644 methods/powell/mnbrak.hpp

(limited to 'methods')

diff --git a/methods/bfgs/bfgs.hpp b/methods/bfgs/bfgs.hpp
new file mode 100644
index 0000000..e5172ef
--- /dev/null
+++ b/methods/bfgs/bfgs.hpp
@@ -0,0 +1,200 @@
+#ifndef BFGS_METHOD
+#define BFGS_METHOD
+#define OPT_HEADER
+#include <core/optimizer.hpp>
+//#include <blitz/array.h>
+#include <limits>
+#include <cstdlib>
+#include <core/opt_traits.hpp>
+#include <math/num_diff.hpp>
+#include <cassert>
+#include <cmath>
+#include <ctime>
+#include <vector>
+#include <algorithm>
+/*
+ *
+*/
+#include <iostream>
+using std::cout;
+using std::endl;
+
+namespace opt_utilities
+{
+
+  template <typename rT,typename pT>
+  class bfgs_method
+    :public opt_method<rT,pT>
+  {
+  private:
+    pT start_point;
+    rT threshold;
+    func_obj<rT,pT>* p_fo;
+    optimizer<rT,pT>* p_optimizer;
+    typedef typename element_type_trait<pT>::element_type element_type;
+    element_type* mem_pool;
+    element_type** invBk;
+    element_type** invBk1;
+    bool bstop;
+  private:
+    const char* do_get_type_name()const
+    {
+      return "asexual genetic algorithm";
+    }
+    
+    rT func(const pT& x)
+    {
+      assert(p_fo!=0);
+      return p_fo->eval(x);
+    }
+
+  public:
+    bfgs_method()
+      :p_fo(0),p_optimizer(0),
+       invBk(0),invBk1(0);
+    {
+      
+    }
+    
+    virtual ~bfgs_method()
+    {     
+      destroy_workspace();
+    };
+    
+    bfgs_method(const bfgs_method<rT,pT>& rhs)
+      :p_fo(rhs.p_fo),p_optimizer(rhs.p_optimizer),
+       threshold(rhs.threshold)
+    {
+    }
+
+    bfgs_method<rT,pT>& operator=(const bfgs_method<rT,pT>& rhs)
+    {
+      p_fo=rhs.p_fo;
+      p_optimizer=rhs.p_optimizer;
+      threshold=rhs.threshold;
+    }
+
+    opt_method<rT,pT>* do_clone()const
+    {
+      return new bfgs_method<rT,pT>(*this);
+    }
+    
+    void init_workspace(int n)
+    {
+      destroy_workspace();
+      mem_pool=new element_type[n*n*2];
+      invBk=new element_type*[n];
+      invBk1=new element_type*[n];
+
+      for(size_t i=0;i!=n;++i)
+	{
+	  invBk[i]=mem_pool+i*n;
+	  invBk1[i]=invBk[i]+n*n;
+	}
+      for(size_t i=0;i!=n;++i)
+	{
+	  for(size_t j=0;j!=n;++j)
+	    {
+	      invBk[i][j]=(i==j?1:0);
+	    }
+	}
+    }
+
+    void destroy_workspace()
+    {
+      delete[] mem_pool;
+      delete[] invBk;
+      delete[] invBk1;
+    }
+
+  public:
+    
+    void do_set_start_point(const pT& p)
+    {
+      start_point=p;
+      init_workspace(get_size(p));
+    }
+
+    pT do_get_start_point()const
+    {
+    }
+    
+    void do_set_precision(rT t)
+    {
+      threshold=t;
+    }
+
+    rT do_get_precision()const
+    {
+      return threshold;
+    }
+
+    void do_set_optimizer(optimizer<rT,pT>& o)
+    {
+      p_optimizer=&o;
+      p_fo=p_optimizer->ptr_func_obj();
+    }
+    
+    pT do_optimize()
+    {
+      pT s;
+      resize(s,get_size(start_point));
+      pT old_grad;
+      pT y;
+      resize(old_grad,get_size(start_point));
+      resize(y,get_size(start_point));
+      for(size_t i=0;i!=get_size(p);++i)
+	{
+	  set_element(old_grad,i,gradient(*p_fo,start_point,i));
+	  set_element(s,i,0);
+	  for(size_t j=0;j!=get_size(p);++j)
+	    {
+	      s[i]+=invBk[i][j]*old_grad[j];
+	    }
+	}
+      linmin(start_point,s,*p_fo);
+      for(size_t i=0;i!=get_size(p);++i)
+	{
+	  set_element(y,gradient(*p_fo,start_point,i)-get_element(old_grad,i));
+	}
+      rT sy=0;
+      pT invBy;
+      pT yinvB;
+      resize(invBy,get_size(p));
+      resize(yinvB,get_size(p));
+      for(size_t i=0;i!=get_size(p);++i)
+	{
+	  sy+=s[i]*y[i];
+	  for(size_t j=0;j!=get_size(p);++j)
+	    {
+	      invBy[i]+=invBk[i][j]*y[j];
+	      yinvB[i]+=y[j]*invBk[j][i];
+	    }
+	}
+      rT yinvBy=0;
+      for(size_t i=0;i!=get_size(p);++i)
+	{
+	  yinvBy+=invBy[i]*y[i];
+	}
+      for(size_t i=0;i<get_size(p);++i)
+	{
+	  for(size_t j=0;j<get_size(p);++j)
+	    {
+	      invBk[i][j]+=((sy+yinvBy)*s[i]*s[j]/(sy*sy)-(invBy[i]*s[j]+s[i]*yinvB[j])/(sy));
+	    }
+	}
+      
+    }
+    
+    void do_stop()
+    {
+      bstop=true;
+    }
+
+  };
+
+}
+
+
+#endif
+//EOF
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
diff --git a/methods/powell/bas_util.hpp b/methods/powell/bas_util.hpp
deleted file mode 100644
index ccdce9e..0000000
--- a/methods/powell/bas_util.hpp
+++ /dev/null
@@ -1,65 +0,0 @@
-#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/powell/brent.hpp b/methods/powell/brent.hpp
deleted file mode 100644
index 1a0ee39..0000000
--- a/methods/powell/brent.hpp
+++ /dev/null
@@ -1,111 +0,0 @@
-#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/powell/linmin.hpp b/methods/powell/linmin.hpp
deleted file mode 100644
index 962b052..0000000
--- a/methods/powell/linmin.hpp
+++ /dev/null
@@ -1,105 +0,0 @@
-#ifndef LINMIN_HPP
-#define LINMIN_HPP
-#define OPT_HEADER
-#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&)
-      :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=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
deleted file mode 100644
index 9871473..0000000
--- a/methods/powell/mnbrak.hpp
+++ /dev/null
@@ -1,82 +0,0 @@
-#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
diff --git a/methods/powell/powell_method.hpp b/methods/powell/powell_method.hpp
index 1e135b0..270c2b2 100644
--- a/methods/powell/powell_method.hpp
+++ b/methods/powell/powell_method.hpp
@@ -10,7 +10,7 @@
 #include <limits>
 #include <cassert>
 #include <cmath>
-#include "linmin.hpp"
+#include "../linmin/linmin.hpp"
 #include <algorithm>
 #include <iostream>
 
-- 
cgit v1.2.2