From 36076e63361d7c3d010f673d267793797de122b8 Mon Sep 17 00:00:00 2001
From: astrojhgu <astrojhgu@ed2142bd-67ad-457f-ba7c-d818d4011675>
Date: Tue, 17 Nov 2009 17:08:23 +0000
Subject: git-svn-id: file:///home/svn/opt_utilities@96
 ed2142bd-67ad-457f-ba7c-d818d4011675

---
 misc/optvec.hpp      | 78 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 statistics/chisq.hpp | 51 +++++++++++++++++++++++++++++++++-
 vmodel/lin1d.hpp     | 58 ++++++++++++++++++++++++++++++++++++++
 3 files changed, 186 insertions(+), 1 deletion(-)
 create mode 100644 vmodel/lin1d.hpp

diff --git a/misc/optvec.hpp b/misc/optvec.hpp
index d23934c..d950c5f 100644
--- a/misc/optvec.hpp
+++ b/misc/optvec.hpp
@@ -1,4 +1,8 @@
+#ifndef OPTVEC_HPP
+#define OPTVEC_HPP
+
 #include <vector>
+#include <cassert>
 
 namespace opt_utilities
 {
@@ -35,4 +39,78 @@ namespace opt_utilities
       return dynamic_cast<const std::vector<T>&>(*this);
     }
   };
+
+  template<typename T>
+  optvec<T> operator+(const optvec<T>& x1,const optvec<T>& x2)
+  {
+    
+    optvec<T> result(x1);
+    for(size_t i=0;i!=x1.size();++i)
+      {
+	result[i]=result[i]+x2.at(i);
+      }
+    return result;
+  }
+
+  template<typename T>
+  optvec<T> operator-(const optvec<T>& x1,const optvec<T>& x2)
+  {
+    
+    optvec<T> result(x1);
+    for(size_t i=0;i!=x1.size();++i)
+      {
+	result[i]=result[i]-x2.at(i);
+      }
+    return result;
+  }
+  
+  template<typename T>
+  optvec<T> operator*(const optvec<T>& x1,const optvec<T>& x2)
+  {
+    
+    optvec<T> result(x1);
+    for(size_t i=0;i!=x1.size();++i)
+      {
+	result[i]=result[i]*x2.at(i);
+      }
+    return result;
+  }
+
+  template<typename T>
+  optvec<T> operator/(const optvec<T>& x1,const optvec<T>& x2)
+  {
+    
+    optvec<T> result(x1);
+    for(size_t i=0;i!=x1.size();++i)
+      {
+	result[i]=result[i]/x2.at(i);
+      }
+    return result;
+  }
+
+  template <typename T>
+  T sum(const optvec<T>& x)
+  {
+    T result=0;
+    for(size_t i=0;i!=x.size();++i)
+      {
+	result+=x[i];
+      }
+    return result;
+  }
+  
+  template<typename T>
+  bool operator<(const optvec<T>& x1,const optvec<T>& x2)
+  {
+    for(size_t i=0;i!=x1.size();++i)
+      {
+	if(x1[i]!=x2[i])
+	  {
+	    return x1[i]<x2[i];
+	  }
+      }
+    return false;
+  }
+  
 };
+#endif
diff --git a/statistics/chisq.hpp b/statistics/chisq.hpp
index 54826d5..f9c2c22 100644
--- a/statistics/chisq.hpp
+++ b/statistics/chisq.hpp
@@ -8,6 +8,7 @@
 #include <core/fitter.hpp>
 #include <iostream>
 #include <vector>
+#include <misc/optvec.hpp>
 #include <cmath>
 using std::cerr;using std::endl;
 
@@ -205,7 +206,55 @@ namespace opt_utilities
       return result;
     }
   };
-#endif  
+#endif
+
+  template<typename T,typename Tp,typename Ts,typename Tstr>
+  class chisq
+    :public statistic<optvec<T>,optvec<T>,Tp,Ts,Tstr>
+  {
+  private:
+    bool verb;
+    int n;
+    
+    typedef optvec<T> Tx;
+    typedef optvec<T> Ty;
+    statistic<Ty,Tx,Tp,Ts,Tstr>* do_clone()const
+    {
+      // return const_cast<statistic<Ty,Tx,Tp>*>(this);
+      return new chisq<Ty,Tx,Tp,Ts,Tstr>(*this);
+    }
+
+    const char* do_get_type_name()const
+    {
+      return "chi^2 statistic";
+    }
+    
+  public:
+    void verbose(bool v)
+    {
+      verb=v;
+    }
+  public:
+    chisq()
+      :verb(false)
+    {}
+    
+    
+
+    Ts do_eval(const Tp& p)
+    {
+      Ts result(0);
+      for(int i=(this->get_data_set()).size()-1;i>=0;--i)
+	{
+	  Ty chi=(this->get_data_set().get_data(i).get_y()-eval_model(this->get_data_set().get_data(i).get_x(),p))/this->get_data_set().get_data(i).get_y_upper_err();
+	  result+=sum(chi*chi);
+
+	}
+      return result;
+    }
+  };
+
+  
 }
 
 #endif
diff --git a/vmodel/lin1d.hpp b/vmodel/lin1d.hpp
new file mode 100644
index 0000000..8ca8074
--- /dev/null
+++ b/vmodel/lin1d.hpp
@@ -0,0 +1,58 @@
+#ifndef VLINEAR_MODEL_H_
+#define VLINEAR_MODEL_H_
+#define OPT_HEADER
+#include <core/fitter.hpp>
+#include <misc/optvec.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+  template <typename T>
+  class lin1d
+    :public model<optvec<T>,optvec<T>,optvec<T>,std::string>
+  {
+    typedef optvec<T> Tv;
+  private:
+    lin1d<T>* do_clone()const
+    {
+      return new lin1d<T>(*this);
+    }
+
+    const char* do_get_type_name()const
+    {
+      return "1d linear model";
+    }
+  public:
+    lin1d()
+    {
+      this->push_param_info(param_info<Tv>("k",1));
+      this->push_param_info(param_info<Tv>("b",0));
+    }
+
+  public:
+    Tv do_eval(const Tv& x,const Tv& param)
+    {
+      Tv result(x.size());
+      
+      //return x*get_element(param,0)+get_element(param,1);
+      for(size_t i=0;i!=x.size();++i)
+	{
+	  result[i]=param[0]*x[i]+param[1];
+	}
+      return result;
+    }
+
+  private:
+    std::string do_get_information()const
+    {
+      return "<math><mrow> <mtext>f(x;k,b)=k x+b</mtext>		\
+    </mrow>								\
+</math>";
+    }
+  };
+}
+
+
+
+#endif
+//EOF
-- 
cgit v1.2.2