aboutsummaryrefslogtreecommitdiffstats
path: root/methods/powell
diff options
context:
space:
mode:
Diffstat (limited to 'methods/powell')
-rw-r--r--methods/powell/bas_util.hpp63
-rw-r--r--methods/powell/brent.hpp110
-rw-r--r--methods/powell/linmin.hpp102
-rw-r--r--methods/powell/mnbrak.hpp81
-rw-r--r--methods/powell/powell_method.hpp251
5 files changed, 607 insertions, 0 deletions
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
Copyright © 2017-2018 Aaron LI. All Rights Reserved.