# include<iostream>// standard input/output # include<vector>// standard vector # include<cppad/cppad.hpp>// the CppAD package
// begin the empty namespace namespace { // define the function Poly(a, x) = a[0] + a[1]*x[1] + ... + a[k-1]*x[k-1] // where x[i] == x^i // "x" and "y" are "Type" type, "a" is a CppAD vector template <classType> Type Poly(const CPPAD_TESTVECTOR(double) &a, const Type &x){ size_t k = a.size(); Type y = 0.; Type x_i = 1.;
for(size_t i = 0; i < k; i++){ y += a[i] * x_i; // y = y + a_i * x^i x_i *= x; // x_i = x_i * x }
return y; } }
intmain(void){
using CppAD::AD; // use AD as abbreviation for CppAD::AD using std::vector; // use vector as abbreviation for std::vector
// 多项式系数 vector of polynomial coefficients size_t k = 5; CPPAD_TESTVECTOR(double) a(k);
for(size_t i = 0; i < k; i++) a[i] = 1.;
// 自变量 domain space vector size_t n = 1; // number of domain space variables vector< AD<double> > ax(n); // vector of domain space variables ax[0] = 3.; // value at which function is recorded
// 因变量 range space vector size_t m = 1; // number of ranges space variables vector< AD<double> > ay(m); // vector of ranges space variables ay[0] = Poly(a, ax[0]); // record operations that compute ay[0]
// 定义 cppAD 的函数映射变量 f // store operation sequence in f: X -> Y and stop recording CppAD::ADFun<double> f(ax, ay);
// compute derivative using operation sequence stored in f vector<double> jac(m * n); // Jacobian of f (m by n matrix) vector<double> x(n); // domain space vector x[0] = 3.; // argument value for computing derivative jac = f.Jacobian(x); // Jacobian for operation sequence
// print the results std::cout << "f'(3) computed by CppAD = " << jac[0] << std::endl;
// check if the derivative is correct int error_code; if( jac[0] == 142. ) error_code = 0; else error_code = 1;
return error_code; }
测试用拟合函数
Gaussian 函数
在本例中，需要拟合的高斯函数的表达式如下：
G(x)=Aπw4ln2e−w24ln2(x−x0)2
其中输入自变量为x，需要估计的参数有(A,x0,w)。
1 2 3 4 5 6 7 8 9 10
template <typename T> T gaussCurve( T A, T x, T x0, T w){ T x2 = (x-x0)*(x-x0); T w2 = w*w; T G = A * (sqrt(4*log(2))/(sqrt(M_PI)*w))*CppAD::exp(-(4*log(2)/w2)*x2); return G; }
// 用 T 应对不同类型的输入和返回，主要是CppAD的类型的输入与返回，用于生成计算图 // 此外，此处必须使用 CppAD::exp() 而不是 std::exp() 以生成计算图
template <typename T> T pseudo_voigt(T A, T n, T x, T x0, T w){ T x2 = (x-x0)*(x-x0); T w2 = w*w; return A*( n*(2.0/M_PI)*(w/(4*x2+w2)) + (1-n)*(sqrt(4*log(2))/(sqrt(M_PI)*w))*CppAD::exp( -(4*log(2)/w2)*x2 ) ); }
template <typename T> voidnls_func(const std::vector<T>& p, std::vector<T>& f, int func_id=0){ // p 是参数向量, f 是残差向量, 包含每个观测点的残差 // 用 T 应对CppAD的类型的输入与返回，用于生成计算图 // 根据传入参数 func_id 确定拟合的函数 T A = p[0]; T x0 = p[1]; T w = p[2];
for (int i = 0; i < f.size(); i++) { T xi = data_x[i]; T yi = data_y[i]; if(func_id){ T n = p[3]; f[i] = yi - pseudo_voigt(A, n, xi, x0, w); }else{ f[i] = yi - gaussCurve(A, xi, x0, w); } } }
voidnls_func(const VectorXd& p, VectorXd& f, int func_id=0){ // 运算符重载，用于实际计算
double A = p[0]; double x0 = p[1]; double w = p[2];
for (int i = 0; i < f.size(); i++) { double xi = data_x[i]; double yi = data_y[i]; if(func_id){ double n = p[3]; f[i] = yi - pseudo_voigt(A, n, xi, x0, w); }else{ f[i] = yi - gaussCurve(A, xi, x0, w); } } }
import pandas as pd import matplotlib.pyplot as plt import math def gaussCurve(A, x, x0, w): x2 = (x - x0) ** 2 w2 = w ** 2 return A * (math.sqrt(4*math.log(2))/(math.sqrt(math.pi)*w)) * math.exp(-(4*math.log(2)/w2)*x2)
def pseudo_voigt(A, n, x, x0, w): x2 = (x - x0) ** 2 w2 = w ** 2 return A * (n*(2.0/math.pi)*(w/(4*x2+w2)) + (1-n)*(math.sqrt(4*math.log(2))/(math.sqrt(math.pi)*w))*math.exp(-(4*math.log(2)/w2)*x2))
def plotCurve(x, y, y_pred, title): plt.scatter(x, y, s=1) plt.plot(x, y_pred, color="r") plt.title(title) plt.show()
if __name__ == "__main__":
# Gaussian 函数 A, x0, w = (100.303, 0.509871, 0.210608) data = pd.read_csv("GaussCurveData.txt", header=None) y_pred = pd.Series([gaussCurve(A,xi,x0,w) for xi in data[0]]) res = y_pred - data[1] title = f"Parameter: $A={A},x_0={x0},\omega={w}$" plotCurve(data[0],data[1],y_pred,title) print("残差平方和",res.pow(2).sum()) print("方差",res.var())
# pseudo-Voigt 函数 A, x0, w, n = (122.06, 0.499887, 0.209804, 0.321467) data = pd.read_csv("pseudo_voigt_data.txt", header=None) y_pred = pd.Series([pseudo_voigt(A,n,xi,x0,w) for xi in data[0]]) res = y_pred - data[1] title = f"Parameter: $A={A},\eta={n},x_0={x0},\omega={w}$" plotCurve(data[0],data[1],y_pred,title) print("残差平方和",res.pow(2).sum()) print("方差",res.var())
