From 3a8b1ca89c9b62e35e24d414d7a79f4c325dfd09 Mon Sep 17 00:00:00 2001
From: astrojhgu <astrojhgu@ed2142bd-67ad-457f-ba7c-d818d4011675>
Date: Sun, 1 Apr 2012 18:17:41 +0000
Subject: git-svn-id: file:///home/svn/opt_utilities@230
 ed2142bd-67ad-457f-ba7c-d818d4011675

---
 core/fp_fo_adapter.hpp                            |  33 ++++
 methods/conjugate_gradient/conjugate_gradient.hpp | 188 ++++++++++++++++++++++
 test/test_cg.cpp                                  |  55 +++++++
 3 files changed, 276 insertions(+)
 create mode 100644 core/fp_fo_adapter.hpp
 create mode 100644 methods/conjugate_gradient/conjugate_gradient.hpp
 create mode 100644 test/test_cg.cpp

diff --git a/core/fp_fo_adapter.hpp b/core/fp_fo_adapter.hpp
new file mode 100644
index 0000000..0c5740b
--- /dev/null
+++ b/core/fp_fo_adapter.hpp
@@ -0,0 +1,33 @@
+#ifndef FP_FO_ADAPTER_HPP
+#define FP_FO_ADAPTER_HPP
+#include "optimizer.hpp"
+
+namespace opt_utilities
+{
+  template <typename rT,typename pT>
+  class fp_fo_adapter
+    :public func_obj<rT,pT>
+  {
+  private:
+    fp_fo_adapter();
+    rT (*fp)(const pT&);
+  public:
+    fp_fo_adapter(rT (*_fp)(const pT&))
+    :fp(_fp)
+    {}
+    
+    fp_fo_adapter<rT,pT>* do_clone()const
+    {
+      return new fp_fo_adapter<rT,pT>(*this);
+    }
+    
+    rT do_eval(const pT& x)
+    {
+      return fp(x);
+    }
+  };
+};
+
+
+#endif
+
diff --git a/methods/conjugate_gradient/conjugate_gradient.hpp b/methods/conjugate_gradient/conjugate_gradient.hpp
new file mode 100644
index 0000000..3fb0f8b
--- /dev/null
+++ b/methods/conjugate_gradient/conjugate_gradient.hpp
@@ -0,0 +1,188 @@
+/**
+   \file conjugate_gradient.hpp
+   \brief powerll optimization method
+   \author Junhua Gu
+ */
+
+#ifndef CONJUGATE_GRADIENT
+#define CONJUGATE_GRADIENT
+#define OPT_HEADER
+#include <core/optimizer.hpp>
+//#include <blitz/array.h>
+#include <limits>
+#include <cassert>
+#include <cmath>
+#include "../linmin/linmin.hpp"
+#include <algorithm>
+#include <iostream>
+#include <math/num_diff.hpp>
+namespace opt_utilities
+{
+  /**
+     \brief Impliment of an optimization method
+     \tparam rT return type of the object function
+     \tparam pT parameter type of the object function
+   */
+  template <typename rT,typename pT>
+  class conjugate_gradient
+    :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;
+    volatile bool bstop;
+    
+    const char* do_get_type_name()const
+    {
+      return "conjugate gradient";
+    }
+  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);
+    }
+
+   
+  private:
+
+
+    
+  public:
+    
+    conjugate_gradient()
+      :threshold(1e-4)
+    {}
+
+    virtual ~conjugate_gradient()
+    {
+    };
+
+    conjugate_gradient(const conjugate_gradient<rT,pT>& rhs)
+      :opt_method<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)
+    {
+    }
+
+    conjugate_gradient<rT,pT>& operator=(const conjugate_gradient<rT,pT>& rhs)
+    {
+      threshold=rhs.threshold;
+      p_fo=rhs.p_fo;
+      p_optimizer=rhs.p_optimizer;
+      start_point=rhs.start_point;
+      end_point=rhs.end_point;
+      threshold=rhs.threshold;
+    }
+    
+    opt_method<rT,pT>* do_clone()const
+    {
+      return new conjugate_gradient<rT,pT>(*this);
+    }
+    
+    void do_set_start_point(const array1d_type& p)
+    {
+      resize(start_point,get_size(p));
+      opt_eq(start_point,p);
+    }
+
+    array1d_type do_get_start_point()const
+    {
+      return start_point;
+    }
+
+    void do_set_lower_limit(const array1d_type& p)
+    {}
+
+    void do_set_upper_limit(const array1d_type& p)
+    {}
+
+    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()
+    {
+      bstop=false;
+      opt_eq(end_point,start_point);
+      pT xn;
+      opt_eq(xn,start_point);
+      pT Delta_Xn1(gradient(*p_fo,start_point));
+      for(size_t i=0;i<get_size(start_point);++i)Delta_Xn1[i]=-Delta_Xn1[i];
+      rT alpha=0;
+      linmin(start_point,Delta_Xn1,alpha,(*p_fo));
+      for(size_t i=0;i<get_size(start_point);++i)xn[i]=start_point[i]+alpha*Delta_Xn1[i];
+      pT LX;
+      opt_eq(LX,Delta_Xn1);
+      for(int n=1;;++n)
+	{
+	  pT Delta_Xn(gradient(*p_fo,xn));
+	  for(size_t i=0;i<get_size(start_point);++i)Delta_Xn[i]=-Delta_Xn[i];
+	  ////calc beta n
+	  rT betan;
+	  rT b1(0),b2(0);
+	  for(size_t i=0;i<get_size(start_point);++i)
+	    {
+	      b1+=Delta_Xn[i]*(Delta_Xn[i]-Delta_Xn1[i]);
+	      b2+=Delta_Xn1[i]*Delta_Xn1[i];
+	    }
+	  betan=max(rT(0),b1/b2);
+	  ////
+	  for(size_t i=0;i<get_size(start_point);++i)
+	    LX[i]=Delta_Xn[i]+betan*LX[i];
+	  linmin(xn,LX,alpha,(*p_fo));
+	  for(size_t i=0;i<get_size(start_point);++i)
+	    xn[i]+=alpha*LX[i];
+	  rT delta=0;
+	  rT xn_abs=0;
+	  for(size_t i=0;i<get_size(start_point);++i)
+	    {
+	      delta+=LX[i]*LX[i];
+	      xn_abs+=xn[i]*xn[i];
+	    }
+	  if(delta*alpha*alpha<threshold)
+	    {
+	      opt_eq(end_point,xn);
+	      break;
+	    }
+	  opt_eq(Delta_Xn1,Delta_Xn);
+	}
+      return end_point;
+    }
+    
+    void do_stop()
+    {
+      bstop=true;
+    }
+
+  };
+
+}
+
+
+#endif
+//EOF
diff --git a/test/test_cg.cpp b/test/test_cg.cpp
new file mode 100644
index 0000000..a080cc5
--- /dev/null
+++ b/test/test_cg.cpp
@@ -0,0 +1,55 @@
+#include <methods/conjugate_gradient/conjugate_gradient.hpp>
+#include <methods/powell/powell_method.hpp>
+#include <vector>
+using namespace std;
+
+using namespace opt_utilities;
+
+class foo
+  :public diff_func_obj<double,vector<double> >
+{
+  foo* do_clone()const
+  {
+    return new foo(*this);
+  }
+
+  double do_eval(const vector<double>& p)
+  {
+    static int call_cnt=0;
+    ++call_cnt;
+    cerr<<call_cnt<<endl;
+    double result=0;
+    for(int i=0;i<p.size();++i)
+      {
+	result+=p[i]*p[i];
+      }
+    return result;
+  }
+
+  vector<double> do_gradient(const vector<double>& p)
+  {
+    static int call_cnt=0;
+    ++call_cnt;
+    cerr<<call_cnt<<endl;
+    vector<double> result(p.size());
+    for(int i=0;i<p.size();++i)
+      {
+	result[i]=2*p[i];
+      }
+    return p;
+  }
+};
+
+
+int main()
+{
+  optimizer<double,vector<double> > opt;
+  //opt.set_opt_method(conjugate_gradient<double,vector<double> >());
+  opt.set_opt_method(powell_method<double,vector<double> >());
+  opt.set_func_obj(foo());
+  std::vector<double> p(2);
+  p[0]=p[1]=4;
+  opt.set_start_point(p);
+  p=opt.optimize();
+  cout<<p[0]<<"\t"<<p[1]<<endl;
+}
-- 
cgit v1.2.2