#include<bits/stdc++.h> 
using namespace std;
vector<double>arr;
vector<vector<double>>train;
vector<vector<double>>test;
double minn=0x3f3f3f,maxn=INT_MIN;
double relu(double x)//激活函数
{
    if(x>0)return x;
    return 0;    
}
class matrix{//定义矩阵类
    private:
        vector<vector<double>>arr;
        int m,n;
    public:
        matrix(){}
        matrix(int _m,int _n):m(_m),n(_n){arr.resize(_m,vector<double>(_n,0));}
        matrix(vector<vector<double>>_a):arr(_a){m=_a.size();n=_a[0].size();}
        vector<vector<double>> get_val(){return arr;}
        void reset() //置零
        {
            for(int i=0;i<m;i++)
            {
                for(int j=0;j<n;j++)
                {
                    arr[i][j]=0;
                }
            }
        }
        void T()//矩阵转置
        {
            vector<vector<double>>t(n,vector<double>(m,0));
            for(int i=0;i<n;i++)
            {
                for(int j=0;j<m;j++)
                {
                    t[i][j]=arr[j][i];
                }
            }
            arr=t;
            swap(n,m);
        }
        void stimulate()//激活函数
        {
            for(int i=0;i<m;i++)
            {
                for(int j=0;j<n;j++)
                {
                    arr[i][j]=relu(arr[i][j]);
                }
            }
        }
        matrix operator+(matrix& _m)//重载操作符
        {
            matrix ans=matrix(m,n);
            for(int i=0;i<m;i++)
            {
                for(int j=0;j<n;j++)
                {
                    ans.arr[i][j]=arr[i][j]+_m.arr[i][j]; 
                }
            }
            return ans;
        }
        matrix operator+(double val)
        {
            matrix ans=matrix(m,n);
            for(int i=0;i<m;i++)
            {
                for(int j=0;j<n;j++)
                {
                    ans.arr[i][j]=arr[i][j]+val;
                }
            }
            return ans;
        }
        matrix operator*(matrix& _m)
        {
            int len1=_m.arr.size(),len2=_m.arr[0].size();
            matrix ans=matrix(m,len2);
            for(int i=0;i<m;i++)
            {
                for(int j=0;j<len2;j++)
                {
                    for(int k=0;k<len1;k++)
                    {
                        ans.arr[i][j]+=_m.arr[k][j]*arr[i][k];
                    }
                }
            }
            return ans;
        }
        matrix operator*(double val)
        {
            matrix ans=matrix(m,n);
            for(int i=0;i<m;i++)
            {
                for(int j=0;j<n;j++)
                {
                    ans.arr[i][j]=arr[i][j]*val;
                }
            }
            return ans;
        }
        matrix operator/(double& val)
        {
            matrix ans=matrix(m,n);
            for(int i=0;i<m;i++)
            {
                for(int j=0;j<n;j++)
                {
                    ans.arr[i][j]=arr[i][j]/val;
                }
            }
            return ans;
        }
        matrix operator-(matrix& _m)
        {
            matrix ans=matrix(m,n);
            for(int i=0;i<m;i++)
            {
                for(int j=0;j<n;j++)
                {
                    ans.arr[i][j]=arr[i][j]-_m.arr[i][j];
                }
            }
            return ans;
        }
        matrix operator-(double& val)
        {
            matrix ans=matrix(m,n);
            for(int i=0;i<m;i++)
            {
                for(int j=0;j<n;j++)
                {
                    ans.arr[i][j]=arr[i][j]-val;
                }
            }
            return ans;
        }
        void display()
        {
            for(int i=0;i<m;i++)
            {
                for(int j=0;j<n;j++)
                {
                    cout<<arr[i][j]<<" ";
                }
                cout<<endl;
            }
        }
};
class NetWork{
    private:
    vector<matrix>w;//权重
    vector<matrix>b;//偏置
    vector<matrix>grand;//权重梯度
    vector<matrix>bgrand;//偏置梯度
    public:
        NetWork(vector<matrix>&_w,vector<matrix>& _b):w(_w),b(_b)
        {
            grand=_w;bgrand=_b;
        }
        
        double forward(int index)//向前计算
        {
            vector<vector<double>>temp(1);
            temp[0]=train[index];temp[0].pop_back();
            matrix m=matrix(temp);
            m=m*w[0];
            m=m+b[0];
            m.stimulate();
            m=m*w[1];
            m=m+b[1];
            m.stimulate();
            return m.get_val()[0][0];
        }
        double loss()//均方差作为损失函数
        {
            double ans=0;
            for(int i=0;i<train.size();i++)
            {
                double f=forward(i)-train[i][13];
                ans+=f*f;
            }
            
            return ans/train.size();
        }
        matrix half_forward(int index)//向前计算一半
        {
            vector<vector<double>>temp(1);
            temp[0]=train[index];temp[0].pop_back();
            matrix m=matrix(temp);
            m=m*w[0];
            m=m+b[0];
            m.stimulate();
            return m;
        }
        void get_grand()//计算梯度
        {
            double len=train.size();
            grand[0].reset();grand[1].reset();
            bgrand[0].reset();bgrand[1].reset();
            for(int i=0;i<train.size();i++)
            {
                double f=forward(i)-train[i][13];//将损失函数对w[1][i][j]求偏导,即可得相应梯度
                f*=2;
                matrix m=half_forward(i);
                m.T();
                m=m*f;
                grand[1]=grand[1]+m;
                bgrand[1]=bgrand[1]+f;
            }
            grand[1]=grand[1]/len;
            bgrand[1]=bgrand[1]/len;        //从后向前更新
            for(int i=0;i<train.size();i++)
            {
                double f=forward(i)-train[i][13];
                vector<vector<double>>temp(1);
                temp[0]=train[i];temp[0].pop_back();
                matrix m=matrix(temp);
                f*=2;
                m=grand[1]*m;
                m.T();
                m=m*f;
                grand[0]=grand[0]+m;
                matrix t=grand[1]*f;
                t.T();
                bgrand[0]=bgrand[0]+t;
            }
            grand[0]=grand[0]/len;
            bgrand[0]=bgrand[0]/len;
        }
        void backward(int times,double eta)//反向传播
        {
            int cnt=0;
            while(cnt<times)
            {
                cnt++;
                get_grand();//更新梯度
                for(int i=0;i<grand.size();i++)//更新权重和偏置
                {
                    matrix m=grand[i]*eta;
                    w[i]=w[i]-m;
                    matrix t=bgrand[i]*eta;
                    b[i]=b[i]-t;
                }
                if(cnt%10==0){cout<<"损失函数："<<loss()<<" "<<endl;}
            }
        }
        double precast(int index)//预测数据
        {
            vector<vector<double>>temp(1);
            temp[0]=test[index];temp[0].pop_back();
            matrix m=matrix(temp);
            m=m*w[0];
            m=m+b[0];
            m.stimulate();
            m=m*w[1];
            m=m+b[1];
            m.stimulate();
            return m.get_val()[0][0];
        }
        vector<matrix> get_w()//返回模型
        {
            return w;
        }
        vector<matrix> get_b()//返回模型
        {
            return b;
        }
};
void normalization(double& data)//归一化范围为0~1
{
    data=(data-minn)/(maxn-minn);
}
void load_data()//读入数据集
{
    ifstream ifp;
    ifp.open("housing.data",ios::in);
    if(!ifp.is_open()){cout<<"operation fails"<<endl;}
    double t;
    while(ifp>>t)//输入数据
    {
        arr.emplace_back(t);minn=min(minn,t);maxn=max(maxn,t);
    }
    ifp.close();
}
void load_mode(vector<matrix>&weigh,vector<matrix>&bias)//保存模型
{
    ifstream ifp;
    ifp.open("mode.txt",ios::in);
    
    for(int i=0;i<weigh.size();i++)
    {
        vector<vector<double>>w;
        w=weigh[i].get_val();
        for(int j=0;j<w.size();j++)
        {
            for(int k=0;k<w[j].size();k++)
            {
                ifp>>w[j][k];
            }
        }
        weigh[i]=matrix(w);
    }
    for(int i=0;i<bias.size();i++)
    {
        vector<vector<double>>b;
        b=bias[i].get_val();
        for(int j=0;j<b.size();j++)
        {
            for(int k=0;k<b[j].size();k++)
            {
                ifp>>b[j][k];
            }
        }
        bias[i]=matrix(b);
    }
    ifp.close();
}
void load_train()//读入训练集
{
    int len=arr.size()/14;
    train.resize(len*8/10);
    int index=0;
    for(int i=0;i<(len*8/10)*14;i++)//一维转二维
    {
        normalization(arr[i]);//数据归一化
        if(i!=0&&(i%14==0))index++;
        train[index].emplace_back(arr[i]);
    }
}
void load_test()//读入测试集
{
    int len=arr.size()/14;
    test.resize(len-(len*8)/10);
    int index=0;
    for(int i=(len*8/10)*14;i<len*14;i++)//一维转二维
    {
        normalization(arr[i]);
        if(i!=(len*8/10)*14&&(i%14==0))index++;
        test[index].emplace_back(arr[i]);
    }
}
void save_mode(vector<matrix>weigh,vector<matrix>bias)//保存模型
{
    ofstream ofp;
    ofp.open("mode.txt",ios::out);
    for(int i=0;i<weigh.size();i++)
    {
        vector<vector<double>>w;
        w=weigh[i].get_val();
        for(int j=0;j<w.size();j++)
        {
            for(int k=0;k<w[j].size();k++)
            {
                ofp<<w[j][k]<<" ";
            }
        }
    }
    for(int i=0;i<bias.size();i++)
    {
        vector<vector<double>>b;
        b=bias[i].get_val();
        for(int j=0;j<b.size();j++)
        {
            for(int k=0;k<b[j].size();k++)
            {
                ofp<<b[j][k]<<" ";
            }
        }
    }
    ofp.close();
}
void init(vector<matrix>&w,vector<matrix>&b)//随机初始化
{
    w.resize(2);b.resize(2);
    srand((unsigned)time(NULL));   
    vector<vector<double>>temp(13,vector<double>(10));
    for(int i=0;i<13;i++)
    {
        for(int j=0;j<10;j++)
        {
            temp[i][j]=(rand()%100)/(double)1000;
        }
    }
    w[0]=matrix(temp);
    vector<vector<double>>t(10,vector<double>(1));
    for(int i=0;i<10;i++)
    {
        for(int j=0;j<1;j++)
        {
            t[i][j]=(rand()%100)/(double)1000;
        }
    }
    w[1]=matrix(t);
    vector<vector<double>>t1(1,vector<double>(10));
    for(int i=0;i<1;i++)
    {
        for(int j=0;j<10;j++)
        {
            t1[i][j]=(rand()%100)/(double)1000;
        }
    }
    b[0]=matrix(t1);
    vector<vector<double>>t2(1,vector<double>(1));
    t2[0][0]=(rand()%100)/(double)1000;
    b[1]=matrix(t2);
}
int main()
{
    vector<matrix>w;//权重
    vector<matrix>b;//偏置
    init(w,b);
    load_data();
    load_train();
    load_test();
    load_mode(w,b);
    NetWork nw=NetWork(w,b);
    nw.backward(100,0.1);//模型训练
    for(int i=0;i<test.size();i++)
    cout<<"实际值:"<<test[i][13]*(maxn-minn)+minn<<"\t"<<"预测值:"<<" "<<nw.precast(i)*(maxn-minn)+minn<<endl;//预测测试数据
    save_mode(nw.get_w(),nw.get_b());
    return 0;
} 

其中文件"mode.txt"和"housing.data"一个是保存参数，一个是开源数据，可以在网上下载