动态规划--分段最小二乘法

#include <iostream>
using namespace std;
#include <fstream>
const int c = 10;
double Min2Method(double **X,int s_num,int e_num);


int main()
{
int n=0,i,j;
cout<<"please input the number of points n = ";
cin>>n;
//initial array arr[i][j]
    double **arr;
arr = new double *[2];
for(i=0;i<n;i++) arr[i] = new double [n];
    ifstream in("test.txt");//打开文件
    //读数据。。
    for( i = 0; i < n; ++i){
        for( j = 0; j < 2; ++j){
double a;
            in >>a;
arr[j][i]=a;
        }
    }
    in.close();//关闭文件
    //输出结果
    for( i = 0; i < n; ++i)
{
        for(int j = 0; j < 2; ++j){
            cout<<arr[j][i]<<" ";
        }
        cout<<endl;
    }
 system("PAUSE");


 //err_ij  definition
double err_ij[100][100];
for( j=0;j<n;j++)
{
for(i=0;i<=j;i++)
{
err_ij[i][j] =Min2Method(arr,i,j);
cout<<"err"<<i<<j<<"="<<err_ij[i][j]<<"\t";
}
cout <<endl;
}
//find the best solution
double *min_err;
min_err = new double[n+1];
min_err[0] = 0;
int *index_err;
index_err = new int [n];
for(i=0;i<n;i++) index_err[i] = 0;
//求最佳解
if( n == 1) {cout<<"please input points more than 1";}
else
{
for(j=1;j<n;j++)
{
double temp = err_ij[0][j]+c;
for(i=1;i<=j;i++)
{
if( err_ij[i][j] +  min_err[i-1] + c < temp )  
{   
index_err[j] = i;
temp = err_ij[i][j] +  min_err[i-1] +c;
}
} 
min_err[j] = temp;
}
}
    
cout<<"**************this is the answer***********************"<<endl<<endl<<endl;
for(i=0;i<n;i++)
{
cout << "min_err"<<i<<"="<<min_err[i]<<" "<<"index_err"<<i<<"="<<index_err[i]+1<<endl;//test
}
    system("PAUSE");
    return 0;
}


//minsquare function
double Min2Method(double **X,int s_num,int e_num)  
{  
    //b_interrupt截距，k_slope斜率，s_num起点，e_num终点，err误差  
    double err=0;
if (s_num==e_num) 
{
err=0;
}
else
{
int i;  
double SumX, SumY, SumXY, SumX2;  
SumX=0;  
SumX2=0;  
for(i=s_num;i<=e_num;i++)  
{  
SumX+=X[0][i];  
SumX2+=(X[0][i]*X[0][i]);  
}
  SumY=0;  
for(i=s_num;i<=e_num;i++)  
{  
SumY+=X[1][i];  
}  
SumXY=0;  
for(i=s_num;i<=e_num;i++)  
{  
SumXY+=(X[0][i]*X[1][i]);  
} 
int nCount=e_num-s_num+1;
double b_interrupt=0,k_slope=0;
b_interrupt=((SumX2*SumY-SumX*SumXY)/(nCount*SumX2-SumX*SumX)); 
k_slope=((nCount*SumXY-SumX*SumY)/(nCount*SumX2-SumX*SumX)); 
//testing slope and b
// cout<<"b="<<b_interrupt;
// cout<<"k="<<k_slope;
//testing
for(i=s_num;i<=e_num;i++) 
{
err+=(X[1][i]-k_slope*X[0][i]-b_interrupt)*(X[1][i]-k_slope*X[0][i]-b_interrupt);
}
}
return err;
}