模仿R语言c++ 向量类c 矩阵类matrix等（持续更新 欢迎指点）

头文件：

class_of_R.h #ifndef CLASS_OF_R_H_#define CLASS_OF_R_H_#include<iostream>using std::cout;using std::endl;using std::cin;using std::istream;using std::ios_base;#include<vector>using std::vector;#include<initializer_list>using std::initializer_list;#include<algorithm>using std::for_each;using std::sort;using std::copy;using std::find;#include<stdexcept>using std::runtime_error;using std::out_of_range;#include<cstdlib>using std::srand;using std::rand;using std::system;#include<ctime>using std::time;#include<cmath>using std::sqrt;using std::log;using std::exp;#include<utility>using std::pair;using std::make_pair;#include<iomanip>using std::setw;#include<set>using std::multiset;using std::set;#include<unordered_map>using std::unordered_map;using std::unordered_multimap;#include<map>using std::multimap;using std::map;#include<string>using std::string;#include<fstream>using std::ifstream;#include<memory>using std::shared_ptr;#include<list>using std::list;class matrix;//类c声明class c{private:    vector<double> c_val;    size_t dim;public:    c():c_val(vector<double>()),dim(0){}    c(initializer_list<double>lst);    c(const vector<double>&c_v):c_val(c_v),dim(c_v.size()){}    c(const matrix&m);    void out(size_t n=1)const{cout<<"The vector is:\n";for(size_t i=0;i<dim;i++)cout<<setw(n)<<c_val[i]<<" ";cout<<endl<<endl;}    double operator[](size_t i)const{return c_val[i];}    double& operator[](size_t i){return c_val[i];}    size_t re_dim()const{return dim;}    friend c operator+(const c&a,const c&b);    friend c operator-(const c&a);    friend c operator-(const c&a,const c&b);    friend c rep(const double&d,size_t n);    friend c operator+(const c&a,const double&d);    friend c operator-(const c&a,const double&d);    friend c operator+(const double&d,const c&a);    friend c operator-(const double&d,const c&a);    friend c operator*(const double&d,const c&a);    friend c operator*(const c&a,const double&d);    friend c operator/(const double&d,const c&a);    friend c operator/(const c&a,const double&d);    friend c operator*(const c&a,const c&b);    friend bool operator==(const c&a,const c&b){return a.c_val==b.c_val;}    friend bool operator!=(const c&a,const c&b){return a.c_val!=b.c_val;}    vector<double> const_re_c_val()const{return c_val;}    //在使用const_re_c_val()时可进行拷贝后使用};//类matrix声明class matrix{private:    vector<c> m_val;    size_t rdim;    size_t cdim;public:    matrix():m_val(vector<c>()),rdim(0),cdim(0){}    matrix(const vector<c>&m):m_val(m),rdim(m_val.size()),cdim(m_val[0].re_dim()){}    matrix(size_t m,size_t n,const double&d):m_val(vector<c>(m,vector<double>(n,d))),rdim(m),cdim(n){}    matrix(const c&v,size_t r,size_t c);    matrix(const c&v);    c operator[](size_t i)const{return m_val[i];}    c& operator[](size_t i){return m_val[i];}    void out(size_t n=3)const;    c col(size_t i)const;    void col_change(size_t i,const c&a);    friend matrix operator+(const matrix&a,const matrix&b);    friend matrix operator-(const matrix&a);    friend matrix operator-(const matrix&a,const matrix&b);    friend matrix operator+(const matrix&a,const double&d);    friend matrix operator+(const double&d,const matrix&a);    friend matrix operator-(const matrix&a,const double&d);    friend matrix operator-(const double&b,const matrix&a);    friend matrix operator*(const matrix&a,const double&d);    friend matrix operator*(const double&d,const matrix&a);    friend matrix operator/(const matrix&a,const double&d);    friend matrix operator/(const double&d,const matrix&a);    void Swap(size_t i,size_t j)    {c temp=this->operator[](i);this->operator[](i)=this->operator[](j);this->operator[](j)=temp;}    void col_Swap(size_t i,size_t j)    {        auto temp=this->col(i);        this->col_change(i,this->col(j));        this->col_change(j,temp);    }    friend bool operator==(const matrix&a,const matrix&b){return a.m_val==b.m_val;}    friend bool operator!=(const matrix&a,const matrix&b){return a.m_val!=b.m_val;}    pair<size_t,size_t> re_dim()const{return make_pair(rdim,cdim);}    friend c operator%(const matrix&m,const c&v);    friend matrix operator%(const matrix&a,const matrix&b);    friend matrix operator%(const c&v,const matrix&m);};//类factor声明class factor{private:    vector<pair<double,size_t>>f_val;    map<double,size_t>f_val_to_levels;    vector<size_t>i_val;    set<double>Levels;public:    factor(const c&v=c());    void out(size_t n=1);    friend c levels(const factor&f);    friend size_t nlevels(const factor&f);    friend factor as_integer(const factor&f);    friend factor as_factor(const c&v);}; #endif // CLASS_OF_R_H_ 类成员函数定义文件：class_of_R.cpp #include"class_of_R.h"//类c的有关定义c::c(initializer_list<double>lst){    vector<double>out;    for_each(lst.begin(),lst.end(),[&](const double&d){out.push_back(d);});    c_val=out;    dim=out.size();}c::c(const matrix&m){    vector<double>temp;        for(size_t j=0;j<m.re_dim().second;j++)                for(size_t i=0;i<m.re_dim().first;i++)        temp.push_back(m[i][j]);    c_val=temp;    dim=temp.size();} //类matrix的有关定义void matrix::out(size_t n)const{    cout<<"The matrix is:\n";    for(size_t i=0;i<rdim;i++)    {        for(size_t j=0;j<cdim;j++)            cout<<setw(n)<<(*this)[i][j]<<" ";        cout<<endl;    }    cout<<endl;}c matrix::col(size_t i)const{    vector<double>temp;    for(size_t i=0;i<rdim;i++)        temp.push_back(0);    c out(temp);    for(size_t j=0;j<rdim;j++)        out[j]=(*this)[j][i];    return out;}void matrix::col_change(size_t i,const c&a){    for(size_t j=0;j<a.re_dim();j++)        (*this)[j][i]=a[j];}matrix::matrix(const c&v,size_t r,size_t c):rdim(r),cdim(c){    matrix out(r,c,0);    for(size_t i=0;i<v.re_dim();i++)        out[i%r][i/r]=v[i];    (this->m_val)=out.m_val;}matrix::matrix(const c&v):rdim(v.re_dim()),cdim(1){    matrix out(v.re_dim(),1,0);    for(size_t i=0;i<out.rdim;i++)        out[i][0]=v[i];    this->m_val=out.m_val;} //类factor有关定义factor::factor(const c&v){    vector<double>temp=v.const_re_c_val();    for_each(temp.begin(),temp.end(),             [&](const double&d){Levels.insert(d);}             );    vector<double>s_to_v(Levels.begin(),Levels.end());    vector<pair<double,size_t>>ttemp;    for(size_t i=0;i<s_to_v.size();i++)        ttemp.push_back(make_pair(s_to_v[i],i));    for(size_t i=0;i<temp.size();i++)    {     for(size_t j=0;j<ttemp.size();j++)     {         if(temp[i]==ttemp[j].first)         {            f_val.push_back(ttemp[j]);         }      }    }    for(size_t i=0;i<ttemp.size();i++)        f_val_to_levels.insert(ttemp[i]);        auto a=f_val.begin();    for(size_t i=0;i<f_val.size();i++)    {         i_val.push_back(a->second);         a++;    }}void factor::out(size_t n){    vector<double>temp(f_val.size(),0);    c out0(temp);    auto a=f_val.begin();    for(size_t i=0;i<out0.re_dim();i++)    {        out0[i]=a->first;        a++;    }    vector<double>temp0(Levels.size(),0);    c out1(temp0);    set<double>::iterator b=Levels.begin();    for(size_t i=0;i<out1.re_dim();i++)    {        out1[i]=*b;        b++;    }    out0.out(n);    cout<<"Levels: ";    cout<<endl;    out1.out(n);}  其他函数即测试程序文件：use_R.cpp #include"class_of_R.h"//若干友元函数c operator+(const c&a,const c&b){    size_t MAX=b.re_dim();    if(a.re_dim()>b.re_dim())        MAX=a.re_dim();    vector<double>temp0(MAX,0);    vector<double>temp1(MAX,0);    vector<double>temp2(MAX,0);    copy(a.c_val.begin(),a.c_val.end(),temp0.begin());    copy(b.c_val.begin(),b.c_val.end(),temp1.begin());    for(size_t i=0;i<MAX;i++)        temp2[i]=temp0[i]+temp1[i];    return c(temp2);}c operator-(const c&a){    c out(a);    for(size_t i=0;i<out.re_dim();i++)        out[i]=-out[i];    return out;}c operator-(const c&a,const c&b){    return (a+(-b));}c rep(const double&d,size_t n){    return c(vector<double>(n,d));}c operator+(const c&a,const double&d){    return a+rep(d,a.re_dim());}c operator-(const c&a,const double&d){    return (a+(-d));}c operator+(const double&d,const c&a){    return a+d;}c operator-(const double&d,const c&a){    return (-a)+d;}c operator*(const double&d,const c&a){    c out(a);    for(size_t i=0;i<a.re_dim();i++)        out[i]*=d;    return out;}c operator*(const c&a,const double&d){    return d*a;}c operator/(const double&d,const c&a){    c out(a);    for(size_t i=0;i<a.re_dim();i++)        out[i]=d/a[i];    return out;}c operator/(const c&a,const double&d){    c out(a);    for(size_t i=0;i<a.re_dim();i++)        out[i]=a[i]/d;    return out;}c operator*(const c&a,const c&b){    if(a.re_dim()!=b.re_dim())        throw runtime_error("dim error!\n");    c out(a);    for(size_t i=0;i<a.re_dim();i++)        out[i]*=b[i];    return out;}matrix operator+(const matrix&a,const matrix&b){    size_t r_MAX=a.rdim,c_MAX=a.cdim;    if(a.rdim<b.rdim)        r_MAX=b.rdim;    if(a.cdim<b.cdim)        c_MAX=b.cdim;    matrix tempa(r_MAX,c_MAX,0);    for(size_t i=0;i<a.rdim;i++)        for(size_t j=0;j<a.cdim;j++)        tempa[i][j]=a[i][j];    matrix tempb(r_MAX,c_MAX,0);    for(size_t i=0;i<b.rdim;i++)        for(size_t j=0;j<b.cdim;j++)        tempb[i][j]=b[i][j];    matrix tempout(r_MAX,c_MAX,0);    for(size_t i=0;i<tempout.rdim;i++)        for(size_t j=0;j<tempout.cdim;j++)        tempout[i][j]=tempa[i][j]+tempb[i][j];    return tempout;}matrix operator-(const matrix&a){    matrix out(a);    for(size_t i=0;i<out.rdim;i++)        for(size_t j=0;j<out.cdim;j++)        out[i][j]=-out[i][j];    return out;}matrix operator-(const matrix&a,const matrix&b){    return ((-b)+a);}matrix operator+(const matrix&a,const double&d){    return a+matrix(a.rdim,a.cdim,d);}matrix operator+(const double&d,const matrix&a){    return a+d;}matrix operator-(const matrix&a,const double&d){    return a+(-d);}matrix operator-(const double&b,const matrix&a){    return (-a)+b;}matrix operator*(const matrix&a,const double&d){    matrix out(a);    for(size_t i=0;i<out.rdim;i++)        for(size_t j=0;j<out.cdim;j++)        out[i][j]*=d;    return out;}matrix operator*(const double&d,const matrix&a){    return a*d;}matrix operator/(const matrix&a,const double&d){    return a*(1/d);}matrix operator/(const double&d,const matrix&a){    matrix out(a.rdim,a.cdim,0);    for(size_t i=0;i<out.rdim;i++)        for(size_t j=0;j<out.cdim;j++)        out[i][j]=d/a[i][j];    return out;}c levels(const factor&f){    vector<double> temp(f.Levels.begin(),f.Levels.end());    c done(temp);    done.out();    return done;}size_t nlevels(const factor&f){    size_t out=f.Levels.size();    cout<<out<<endl;    return out;}//返回输入因子相应在存储中的下标构成的因子factor as_integer(const factor&f){    vector<double> temp(f.i_val.begin(),f.i_val.end());    factor temp0=factor(c(temp));    for_each(temp0.f_val.begin(),temp0.f_val.end(),             [](const pair<double,size_t>&p)             {                 cout<<p.first<<" ";             }             );             cout<<endl;    return temp0;}factor as_factor(const c&v){    factor temp0=factor(v);    temp0.out();    return temp0;}//基本运算函数matrix diag(const c&v){    matrix out(v.re_dim(),v.re_dim(),0);    for(size_t i=0;i<v.re_dim();i++)        out[i][i]=v[i];    return out;}matrix diag(const matrix&m){    if(m.re_dim().first!=m.re_dim().second)        throw runtime_error("dim error!\n");    matrix out(m.re_dim().first,m.re_dim().second,0);    for(size_t i=0;i<out.re_dim().first;i++)        out[i][i]=m[i][i];    return out;}double inter_product(const c&a,const c&b){    if(a.re_dim()!=b.re_dim())        throw runtime_error("dim error!\n");    double out=0;    for(size_t i=0;i<a.re_dim();i++)    out+=(a[i]*b[i]);    return out;}c operator%(const matrix&m,const c&v){    if(m.re_dim().second!=v.re_dim())        throw runtime_error("dim error!\n");    c out(rep(0,m.re_dim().first));    for(size_t i=0;i<out.re_dim();i++)        out[i]=inter_product(m[i],v);    return out;}//同下double prod(const c&v){    double out=1;    for(size_t i=0;i<v.re_dim();i++)        out*=v[i];    return out;}//对于matrix的求和由c复制构造函数的重载函数定义：c(const matrix&m)double sum(const c&v){    double out=v[0];    for(size_t i=1;i<v.re_dim();i++)        out+=v[i];    return out;}matrix operator%(const c&v,const matrix&m){    if(v.re_dim()!=m.re_dim().first)        throw runtime_error("dim error!\n");    matrix out(1,m.re_dim().second,0);    for(size_t i=0;i<out.re_dim().second;i++)        out[0][i]=inter_product(v,m.col(i));    return out;}matrix outer_product(const c&a,const c&b){    matrix out(a.re_dim(),b.re_dim(),0);    for(size_t i=0;i<out.re_dim().first;i++)        for(size_t j=0;j<out.re_dim().second;j++)        out[i][j]=a[i]*b[j];    return out;}size_t length(const c&v){    return v.re_dim();}c sort(const c&v){    vector<double>out;    for(size_t i=0;i<v.re_dim();i++)        out.push_back(v[i]);    sort(out.begin(),out.end());    return c(out);}double mean(const c&v){    return sum(v)/length(v);}//其转换版本返回矩阵转为向量的方差double var(const c&v){    return mean((v-mean(v))*(v-mean(v)));}c connect(initializer_list<c> lst){    vector<double>temp;    for_each(lst.begin(),lst.end(),[&](const c&v){for(auto a:v.const_re_c_val())temp.push_back(a);});    return c(temp);}//标准差template<typename any>double STD(const any&a){    return sqrt(var(a)*(length(a))/(length(a)-1));}matrix rbind(initializer_list<c>lst){    size_t r=0;    auto a=lst.begin();    while(a!=lst.end())    {        a++;        r++;    }    size_t co=0;    for_each(lst.begin(),lst.end(),[&](const c&v){if(v.re_dim()>co)co=v.re_dim();});    matrix out(r,co,0);    size_t i=0;    a=lst.begin();    while(a!=lst.end())    {        vector<double>temp(co,0);        for(size_t i=0;i<a->const_re_c_val().size();i++)            temp[i]=a->operator[](i);        out[i]=c(temp);        i++;        a++;    }    return out;}size_t dim(const c&v){    return v.re_dim();}pair<size_t,size_t> dim(const matrix&m){    return m.re_dim();}double min(const c&v){    double temp=v[0];    for(size_t i=1;i<length(v);i++)        if(v[i]<temp)        temp=v[i];    return temp;}double max(const c&v){    double temp=v[0];    for(size_t i=1;i<length(v);i++)        if(v[i]>temp)        temp=v[i];    return temp;}c pmin(initializer_list<c>lst){    c temp=rep(0,lst.size());    auto a=lst.begin();    for(size_t i=0;i<temp.re_dim();i++)    {        temp[i]=min(*a);        a++;    }    return temp;}c pmax(initializer_list<c>lst){    c temp=rep(0,lst.size());    auto a=lst.begin();    for(size_t i=0;i<temp.re_dim();i++)    {        temp[i]=max(*a);        a++;    }    return temp;}template<typename any>pair<double,double> range(const any&v){    cout<<"MIN: "<<min(v)<<" MAX: "<<max(v)<<endl;    return make_pair(min(v),max(v));}double median(const c&v){    return sort(v)[length(v)/2];}double var(const c&a,const c&b){    if(length(a)!=length(b))        throw runtime_error("dim error!\n");    c t1(a-mean(a)),t2(b-mean(b));    return inter_product(t1,t2)/(length(a)-1);}double cov(const c&a,const c&b){    return var(a,b);}double cor(const c&a,const c&b){    return var(a,b)/STD(a)/STD(b);}matrix var(const matrix&a,const matrix&b){    if(a.re_dim().first!=b.re_dim().first)        throw runtime_error("dim error!\n");    matrix temp(a.re_dim().second,b.re_dim().second,0);    for(size_t i=0;i<temp.re_dim().first;i++)        for(size_t j=0;j<temp.re_dim().second;j++)        temp[i][j]=var(a.col(i),b.col(j));    return temp;}matrix cov(const matrix&a,const matrix&b){    return var(a,b);}matrix cor(const matrix&a,const matrix&b){    matrix temp=var(a,b);    for(size_t i=0;i<temp.re_dim().first;i++)        for(size_t j=0;j<temp.re_dim().second;j++)        temp[i][j]/=((STD(a.col(i))*STD(b.col(j))));    return temp;}matrix cov(const matrix&a){    return var(a,a);}matrix cor(const matrix&a){    return cor(a,a);}//由于关联容器的删除方法利用误差排序的方法//利用顺序容器迭代器的删除方法有困难c rank(const c&v){    vector<double>temp00(v.const_re_c_val());    for(size_t i=0;i<temp00.size();i++)        temp00[i]+=(0.00000001*i);    vector<double>temp0(temp00);    multiset<double>temp2(temp0.begin(),temp0.end());    vector<double>temp1(v.re_dim());    size_t i=0,j=0;    double pre_val=min(c(vector<double>(temp2.begin(),temp2.end())));    vector<double>::iterator s0;    while(i<v.re_dim())    {    double val=min(c(vector<double>(temp2.begin(),temp2.end())));    s0=find(temp0.begin(),temp0.end(),val);    if(sqrt((pre_val-val)*(pre_val-val))>0.000001)        j++;    temp1[s0-temp0.begin()]=j;    temp2.erase(val);    i++;    pre_val=val;    }    return  c(temp1);}c rev(const c&v){    vector<double> temp;    auto a=v.const_re_c_val().rbegin();    while(a!=v.const_re_c_val().rend())    {        temp.push_back(*a);        a++;    }    return c(temp);}//向量先排序再输出排序前的下标值c order(const c&v){    unordered_multimap<double,double>m0;    for(size_t i=0;i<length(v);i++)        m0.insert(make_pair(v[i],i));    multimap<double,double>m1(m0.begin(),m0.end());    c temp(rep(0,length(v)));    auto a=m1.begin();    size_t i=0;    while(a!=m1.end())    {        temp[i]=a->second;        i++;        a++;    }    return temp;}c cumsum(const c&v){    c temp(v);    for(size_t i=1;i<dim(v);i++)        temp[i]+=temp[i-1];    return temp;}c cumprod(const c&v){    c temp(v);    for(size_t i=1;i<length(v);i++)        temp[i]*=temp[i-1];    return temp;}matrix cumsum(const matrix&m){    c temp(m);    for(size_t i=1;i<length(m);i++)        temp[i]+=temp[i-1];    return matrix(temp,m.re_dim().first,m.re_dim().second);}matrix cumprod(const matrix&m){    c temp(m);    for(size_t i=1;i<length(m);i++)        temp[i]*=temp[i-1];    return matrix(temp,m.re_dim().first,m.re_dim().second);} //矩阵运算函数//初等行变换化上三角求解行列式double det(const matrix&m){auto f=[](const matrix&mm)->pair<matrix,double>{    if(mm.re_dim().first!=mm.re_dim().second)       throw runtime_error("dim error!\n");    size_t MIN=mm.re_dim().first;    double mult=1;    size_t col=0;    size_t bar=0;    size_t b_c=0;    matrix t(mm);    while(b_c<MIN)    {    for(size_t i=bar+1;i<MIN;)    {        if(t[i][col]==0)            i++;        else        {           t.Swap(i,bar);           mult*=(-1);           break;        }    }    if(bar!=MIN+1)    {    if(t[bar][col]!=0)    {       for(size_t i=bar+1;i<MIN;i++)     {        t[bar][i]=t[bar][i]/t[bar][col];        if(i==bar+1)        mult*=t[bar][col];     }        if(col!=MIN-1)        t[bar][col]=1;    }     for(size_t i=bar+1;i<MIN;i++)     {         t[i]=t[i]-t[i][col]*t[bar];     }    }      col++;      bar++;      b_c=bar;      if(col>bar);      b_c=col;    }return make_pair(t,mult);};    pair<matrix,double> tt=f(m);    return tt.second*prod(diag(tt.first)%rep(1,m.re_dim().first));}//Cramer法则解线性方程组：要求自变量矩阵为方阵且相应行列式不为0c solve(const matrix&m,const c&v){    c temp(rep(0,v.re_dim()));    double low=det(m);    if(m.re_dim().first!=m.re_dim().second)      throw runtime_error("dim error!\n");    if(low==0)      throw runtime_error("singular!\n");    for(size_t i=0;i<v.re_dim();i++)    {        matrix temp0(m);        temp0.col_change(i,v);        temp[i]=det(temp0)/low;    }    return temp;}double tr(const matrix&m){    if(m.re_dim().first!=m.re_dim().second)        throw runtime_error("dim error!\n");    return inter_product(diag(m)%rep(1,m.re_dim().first),rep(1,m.re_dim().first));}matrix T(const matrix&m){    matrix out(m.re_dim().second,m.re_dim().first,0);    for(size_t i=0;i<out.re_dim().first;i++)        for(size_t j=0;j<out.re_dim().second;j++)        out[i][j]=m[j][i];    return out;}matrix cbind(initializer_list<c>lst){    return T(rbind(lst));}matrix operator%(const matrix&a,const matrix&b){    matrix out(a.re_dim().first,b.re_dim().second,0);    for(size_t i=0;i<out.re_dim().first;i++)        for(size_t j=0;j<out.re_dim().second;j++)        out[i][j]=inter_product(a[i],b.col(j));    return out;}double F_norm(const matrix&m){    return sqrt(tr(T(m)%m));} //a,b为区间 d为步长 e为精度 默认为求整数精度：以1为精度特征值 如所求特征值太少可通过减少d改善//只能求实值//可以用F_norm作为a、b的界（-F_norm(m),F_norm(m)）可确保值在范围中//可用于求解但效率低c eigen(const matrix&m,const double&a,const double &b,const double &d=0.0001,const double &e=1){    if(m.re_dim().first!=m.re_dim().second)        throw runtime_error("dim error!\n");    vector<double>temp;    double pre_beg=100000000;    for(double beg=a;beg<b;beg+=d)    {        matrix mm0(diag(rep(beg,m.re_dim().first))-m);        double s0=det(mm0);        if(sqrt(s0*s0)<e)        {            if(sqrt((pre_beg-beg)*(pre_beg-beg))>1)            {             temp.push_back(beg);             pre_beg=beg;             cout<<"Done!\n";             cout<<beg<<endl;//eigenvalue            }        }    }    return c(temp);}//重载求eign的函数 d:步长 可保证绝对求出所有特征值（利用循环增大精度）c eigen(const matrix&m,const double&d=0.1){    double dd=d;    double ra=F_norm(m);    double trr=tr(m);    c temp;    double r=0;    do    {        dd*=0.1;//原来设为0.1        temp=eigen(m,-ra,ra,dd);        if(temp.re_dim()==m.re_dim().first)            break;        r=sum(temp)-trr;    }    while(sqrt((r*r)>1));        //与迹的差距（精度）设定为1 绝对值小于此值的特征值被忽略        //此种方法的一个缺陷为异号特征值的抵消 对于同号特征值有把握        //只适用于所有特征值为实值的情况        //如有复特征值，只输出实数的并循环不能结束    return temp;}//矩阵的余子阵matrix M(const matrix&m,size_t k,size_t l){    matrix temp(m);    set<size_t> s1;    for(size_t i=0;i<m.re_dim().first;i++)        s1.insert(i);    set<size_t> s2;    for(size_t i=0;i<m.re_dim().second;i++)        s2.insert(i);    s1.erase(k);    s2.erase(l);    vector<double>temp00;    for(auto j:s2)        for(auto i:s1)        temp00.push_back(temp[i][j]);    return matrix (c(temp00),m.re_dim().first-1,m.re_dim().second-1);}//方阵的伴随矩阵matrix A_accompany(const matrix&m){    matrix out(m.re_dim().first,m.re_dim().second,0);    for(size_t i=0;i<out.re_dim().first;i++)        for(size_t j=0;j<out.re_dim().second;j++)        out[i][j]=det(M(m,i,j))*(((i+j)%2)?-1:1);        //代数余子式    return T(out);}//求方阵的逆矩阵（利用伴随矩阵方法）matrix solve(const matrix&m){    double de=det(m);    if(sqrt(de*de)<0.000001)     throw runtime_error("singular!\n");    return A_accompany(m)/de;}//利用逆矩阵法求解方程组 要求自变量阵为方阵 且可逆c solve_eq(const matrix&m,const c&v){    c out;    double de=det(m);    if(sqrt(de*de)<0.000001)      throw runtime_error("singular!\n");    return solve(m)%v;}//化列向量组构成的矩阵为正交化后列向量组构成的矩阵（Schmidt 正交化）matrix orth_normal(const matrix &m){    matrix out(rep(0,m.re_dim().first*m.re_dim().second),m.re_dim().first,m.re_dim().second);    vector<c>dec;    dec.push_back(m.col(0)/sqrt(inter_product(m.col(0),m.col(0))));    for(size_t i=1;i<m.re_dim().second;i++)    {        c temp=m.col(i);        for(size_t j=1;j<=i;j++)        temp=temp-inter_product(dec[j-1],m.col(i))*dec[j-1];        temp=temp/sqrt(inter_product(temp,temp));        dec.push_back(temp);    }    for(size_t i=0;i<m.re_dim().second;i++)        out.col_change(i,dec[i]);    return out;}//生成相应位数的由0/1字符构成的字符串vector<string> generator0_1(size_t n){    list<string>out;    list<string>require;    out.push_back("0");    out.push_back("1");    size_t i=0;    while(i<n-1)    {        if(i>0)         out=require;         size_t s0=out.size();    list<string>::iterator pre_begin=out.begin();    auto pre_end=out.end();    size_t k=0;    while(pre_begin!=pre_end)    {        pre_begin++;        k++;    }    pre_begin=out.begin();         for(size_t i=0;i<k;i++)      {          string temp0=*pre_begin+"0";          string temp1=*pre_begin+"1";         out.push_back(temp0);         out.push_back(temp1);         pre_begin++;      }      require.clear();      list<string>::iterator a=out.begin();      for(size_t i=0;i<s0;i++)        a++;      while(a!=out.end())      {          require.push_back(*a);          a++;      }i++;    }   return vector<string>(require.begin(),require.end());}//利用矩阵向量方法生成上述数组matrix matrix_generator0_1(size_t n){    matrix temp0(2,1,0);    matrix temp1;    temp0[0][0]=0;    temp0[1][0]=1;    size_t i=0;    while(i<n-1)    {        if(i>0)         temp0=temp1;//         size_t s0=temp0.re_dim().second;    size_t pre_begin=0;    auto pre_end=temp0.re_dim().first;    temp1=matrix(temp0.re_dim().first*2,temp0.re_dim().second+1,0);    size_t temp1_index=0;        while(pre_begin!=pre_end)      {         vector<double> tt0=temp0[pre_begin].const_re_c_val();         vector<double> tt1=tt0;         tt0.push_back(0);         tt1.push_back(1);         temp1[temp1_index]=c(tt0);         temp1_index++;         temp1[temp1_index]=c(tt1);         temp1_index++;         pre_begin++;      }      i++;    }   return temp1;}//将固定位数的下标值与0/1数值构成的map以vector的形式返回//每个map都是一种组合vector<map<size_t,size_t>> index_map_generator(size_t n){    vector<string>s0(generator0_1(n));    vector<map<size_t,size_t>>outer;    for(size_t i=0;i<s0.size();i++)    {        map<size_t,size_t>tp;        for(size_t j=0;j<s0[i].size();j++)       {        size_t temp=1;        if(s0[i][j]=='0')        temp=0;        tp.insert(make_pair(j,temp));       }        outer.push_back(tp);    }    return outer;}//上述产生方式的矩阵版本vector<map<size_t,size_t>> index_map_matrix_generator(size_t n){    matrix temp=matrix_generator0_1(n);    vector<map<size_t,size_t>> out(temp.re_dim().first);    for(size_t i=0;i<out.size();i++)    {        map<size_t,size_t>m;        for(size_t j=0;j<temp.re_dim().second;j++)            m.insert(make_pair(j,temp[i][j]));        out[i]=m;    }    return out;}//返回上述index_map_generator(m)中指向0_1和为n的子集vector<map<size_t,size_t>> require_dim(size_t m,size_t n){    vector<map<size_t,size_t>>v0=index_map_generator(m);    vector<map<size_t,size_t>>out;    auto a=v0.begin();    while(a!=v0.end())    {        size_t temp=0;        for_each(a->begin(),a->end(),[&](const pair<size_t,size_t>&p){                 if(p.second==1)                    temp++;                 });        if(temp==n)            out.push_back(*a);        a++;    }    return out;}//matrix 所有n阶子式构成的向量vector<double> seq_of_subdet(const matrix&m,size_t n){    vector<double>out;    vector<map<size_t,size_t>>r_map=require_dim(m.re_dim().first,n);    vector<map<size_t,size_t>>c_map=require_dim(m.re_dim().second,n);    for(size_t i=0;i<r_map.size();i++)        for(size_t j=0;j<c_map.size();j++)    {        set<size_t>r_require;        for_each(r_map[i].begin(),r_map[i].end(),[&](const pair<size_t,size_t>&p){if(p.second==1)r_require.insert(p.first);});        set<size_t>c_require;        for_each(c_map[j].begin(),c_map[j].end(),[&](const pair<size_t,size_t>&p){if(p.second==1)c_require.insert(p.first);});        matrix r_matrix(n,m.re_dim().second,0);        auto a=r_require.begin();        for(size_t i=0;i<r_matrix.re_dim().first;i++)        {            r_matrix[i]=m[*a];            a++;        }        matrix c_matrix(r_matrix.re_dim().first,n,0);        auto b=c_require.begin();        for(size_t i=0;i<n;i++)        {            c_matrix.col_change(i,r_matrix.col(*b));            b++;        }        out.push_back(det(c_matrix));    }    return out;}//利用子式的方法求出矩阵的秩size_t rank(const matrix&m){    size_t MAX=m.re_dim().first;    if(m.re_dim().second<m.re_dim().first)        MAX=m.re_dim().second;    size_t out=MAX;    size_t pre_num_of_non0=0;    for(size_t i=MAX;i>0;)    {        vector<double>temp=seq_of_subdet(m,i);        auto a=temp.begin();        while(a!=temp.end())        {            *a=sqrt((*a)*(*a));            a++;        }        size_t num_of_non0=0;        for_each(temp.begin(),temp.end(),[&](const double&d){if(d>0.00000001)num_of_non0++;});        if(pre_num_of_non0==0)            i--;        if(num_of_non0!=0)            return i+1;        pre_num_of_non0=num_of_non0;    }    return out;}//回归系数的求解 要求增广后的数据阵（含常系数）可逆c lm(const c&y,const matrix&m){    matrix X=matrix(m.re_dim().first,m.re_dim().second+1,0);    X.col_change(0,rep(1,X.re_dim().first));    for(size_t i=0;i<m.re_dim().second;i++)        X.col_change(i+1,m.col(i));    if(det(T(X)%X)<0.00000001)        throw runtime_error("singular!\n");    c out(solve(T(X)%X)%T(X)%y);    cout<<"The intercept is: "<<out[0]<<endl;    cout<<"Others are:\n";    for(size_t i=1;i<out.re_dim();i++)        cout<<"x"<<i<<" is "<<out[i]<<endl;    cout<<endl;    return out;}//抽样函数 参数为向量v 抽取的个数n 是否有放回replace 抽取的概率prob//使用此方法前必须在函数调用外设定随机数种子//要求没有0概率点c sample(const c&v,size_t n,bool replace=1,const c&prob=c()){    if(!replace&&v.re_dim()<n)        throw out_of_range("n is too large\n");        //抽样概率    vector<double>probb;    if(prob==c())        probb=rep(1/(double)v.re_dim(),v.re_dim()).const_re_c_val();    else    {        probb=prob.const_re_c_val();        double total=0;        for(auto a:probb)            total+=a;        for(auto &a:probb)            a/=total;    }    ifstream fin("h://runif0.txt");    vector<double>runif;    double val;    while(fin>>val)    {        runif.push_back(val);    }    fin.clear();    //c(probb).out();    c cprobb=cumsum(c(probb));    //cprobb.out();    c require(rep(0,n));    if(replace)    {          size_t i=0;          while(i<n)          {          double uniform=runif[rand()%30000];          //cout<<uniform<<endl;          double conclusion=v[v.re_dim()-1];          for(size_t i=0;cprobb[i]<1;i++)          {              if(uniform<cprobb[i])              {                conclusion=v[i];                break;              }          }          require[i]=conclusion;          i++;          }    }    else        {            size_t i=0;            map<double,double>mdd;            for(size_t i=0;i<v.re_dim();i++)                mdd.insert(make_pair(v[i],c(probb)[i]));            while(i<n)            {            vector<double> temp_v;            vector<double> temp_prob;            for_each(mdd.begin(),mdd.end(),[&](const pair<double,double>&p){                    temp_v.push_back(p.first);                    temp_prob.push_back(p.second);                     });            double d=sample(c(temp_v),1,1,c(temp_prob)).const_re_c_val()[0];            require[i]=d;            mdd.erase(d);            i++;            }        }return require;}//产生几何分布随机数 n 产生个数 p成功概率c rgeom(const size_t& n,const double& p){    vector<double>out;    size_t i=0;    srand(0);    ifstream fin("h://runif0.txt");    vector<double>runif;    double val;    while(fin>>val)    {        runif.push_back(val);    }    fin.clear();    while(i<n)    {        double temp=runif[rand()%30000];        out.push_back(size_t(log(temp)/log(1-p))+1);        i++;    }    return c(out);}//possion分布随机数 n 生成个数 d 分布参数c rpois(const size_t&n,const double&d){    size_t i=0;    vector<double>conclusion;    srand(time(0));    ifstream fin("h://runif0.txt");    vector<double>runif;//储存30000个均匀随机数的vector    double val;    while(fin>>val)    {        runif.push_back(val);    }    fin.clear();    while(i<n)    {        double p=exp(-1*d);        size_t out=0;        double F=p;        double temp=runif[rand()%30000];         while(temp>=F)         {            p=d/(double)(out+1)*p;            F+=p;            out++;         }        conclusion.emplace_back(out);        i++;    }    //for_each(conclusion.begin(),conclusion.end(),[](const double&s){cout<<s<<" ";});    //cout<<endl;    return c(conclusion);}//利用间距生成向量c seq(const double&from,const double&to,const double&by){    if(from>to)        throw runtime_error("from > to\n");    vector<double>out;    auto begin=from;    while(begin<=to)    {        out.push_back(begin);        begin+=by;    }    return c(out);}//生成列表向量元素的所有组合并将结果以矩阵形式输出//前面的generate_grid_ 函数为铺垫而设 expand_grid函数接受c为参数给出输出matrix generate_grid_0(pair<double,c>p0){    matrix out(p0.second.re_dim(),2,0);    for(size_t i=0;i<out.re_dim().first;i++)    {        out[i][0]=p0.first;        out[i][1]=p0.second[i];    }    return out;}matrix generate_grid_1(const c&c0,const c&c1){    vector<double> temp0=c0.const_re_c_val();    size_t len=temp0.size();    matrix out(len*c1.re_dim(),2,0);    for(size_t i=0;i<len;i++)    {        matrix temp00=generate_grid_0(make_pair(temp0[i],c1));       for(size_t j=0;j<c1.re_dim();j++)       {         out[i*c1.re_dim()+j]=temp00[j];       }    }    return out;}matrix generate_grid_2(const matrix&m,const c&v){    matrix temp=generate_grid_1(m.col(m.re_dim().second-1),v);    //temp.out();    //system("pause");    matrix out(temp.re_dim().first,m.re_dim().second+1,0);    //out.out();    //system("pause");    matrix add(m.re_dim().first,m.re_dim().second-1,0);    for(size_t i=0;i<add.re_dim().second;i++)        add.col_change(i,m.col(i));    //add.out();    //system("pause");    for(size_t i=0;i<m.re_dim().first;i++)        for(size_t j=0;j<v.re_dim();j++)        {            for(size_t k=0;k<add.re_dim().second;k++)             out[i*v.re_dim()+j][k]=add[i][k];        }    out.col_change(out.re_dim().second-2,temp.col(0));    out.col_change(out.re_dim().second-1,temp.col(1));    return out;}//若干向量其元素的所有组合matrix expand_grid(initializer_list<c>i0){    vector<c>temp(i0);    if(temp.size()==1)        return matrix(temp[0]);    matrix beg0=generate_grid_1(temp[0],temp[1]);    if(temp.size()==2)        return beg0;    for(size_t i=2;i<temp.size();i++)    {        //beg0.out();        beg0=generate_grid_2(beg0,temp[i]);        //cout<<"Done!\n";        //beg0.out();        //cout<<"Done!\n";    }    return beg0;}//得到向量v的所有n个元素的组合 输出以矩阵排列 每列为一个组合//调用了求得下标组合的 require_dim()//当组合项数太多时主张利用转置按行查看matrix combn(const c&v,size_t n){   auto a=require_dim(v.re_dim(),n);   auto b=vector<set<size_t>>(a.size());   auto a0=a.begin();   auto b0=b.begin();   while(a0!=a.end())   {       for_each(a0->begin(),a0->end(),                [&](const pair<size_t,size_t>&p)                {                    if(p.second==1)                    b0->insert(p.first);                });                a0++;                b0++;   }    auto cc=vector<vector<double>>(b.size());    for(size_t i=0;i<cc.size();i++)    {        cc[i]=vector<double>(b[i].begin(),b[i].end());    }      matrix out(cc[0].size(),cc.size(),0);      for(size_t i=0;i<out.re_dim().second;i++)      {         out.col_change(i,cc[i]);      }      matrix oout(out.re_dim().first,out.re_dim().second,0);      for(size_t i=0;i<oout.re_dim().first;i++)        for(size_t j=0;j<oout.re_dim().second;j++)        oout[i][j]=v[out[i][j]];      return oout;}//上述函数的重载版本 返回0:m的相应n组合matrix combn(size_t m,size_t n){    vector<double>temp;    for(size_t i=0;i<m+1;i++)        temp.push_back(i);    return combn(temp,n);}//向量相应元素在前面是否出现//出现显示1c duplicated(const c&v){    //函数a返回向量cc的前n个元素构成的向量    auto a=[](const c&cc,size_t n)->c    {        c out(rep(0,n));        for(size_t i=0;i<n;i++)            out[i]=cc[i];        return out;    };    vector<c>vc0(v.re_dim());    for(size_t i=0;i<vc0.size();i++)        vc0[i]=a(v,i+1);    vector<vector<double>>vs0(vc0.size());    for(size_t i=0;i<vs0.size();i++)        vs0[i]=vector<double>(vc0[i].const_re_c_val());    c out(rep(0,v.re_dim()));    for(size_t i=1;i<v.re_dim();i++)    {         if(find(vs0[i-1].begin(),vs0[i-1].end(),v[i])!=vs0[i-1].end())         out[i]=1;    }    return out;}//实现unique的操作看似与factor的levels相同//实际在序上是不同的//levels保持double的序（由set有序决定）//unique应保持原向量的序//unique由duplicated定义c unique(const c&v){    c d0=duplicated(v);    vector<pair<double,double>>vm0;    for(size_t i=0;i<v.re_dim();i++)        vm0.push_back(make_pair(v[i],d0[i]));    vector<double>out;    for(size_t i=0;i<vm0.size();i++)    {        if(vm0[i].second==0)            out.push_back(vm0[i].first);    }    return c(out);}//矩阵叉乘 T(a)%b (向量内积的矩阵版本)matrix crossprod(const matrix&a,const matrix&b){    return T(a)%b;}//outer 生成外积 相应的可将外机作用于函数matrix outer(const c&v1,const c&v2,function<double(const double&,const double&)>p=function<double(const double&,const double&)>()){    if(!p)    {        matrix out(v1.re_dim(),v2.re_dim(),0);        for(size_t i=0;i<out.re_dim().first;i++)            for(size_t j=0;j<out.re_dim().second;j++)            out[i][j]=v1[i]*v2[j];        return out;    }    else    {        matrix out(v1.re_dim(),v2.re_dim(),0);        for(size_t i=0;i<out.re_dim().first;i++)            for(size_t j=0;j<out.re_dim().second;j++)            out[i][j]=p(v1[i],v2[j]);        return out;    }}//矩阵Kronecker积matrix kronecker(const matrix&a,const matrix&b){    matrix out(a.re_dim().first*b.re_dim().first,a.re_dim().second*b.re_dim().second,0);    vector<vector<matrix>>vvm0(a.re_dim().first,vector<matrix>(a.re_dim().second));    for(size_t i=0;i<a.re_dim().first;i++)        for(size_t j=0;j<a.re_dim().second;j++)    {       vvm0[i][j]=a[i][j]*b;    }    for(size_t i=0;i<vvm0.size();i++)        for(size_t j=0;j<vvm0[0].size();j++)    {        matrix temp=vvm0[i][j];        for(size_t k=0;k<temp.re_dim().first;k++)            for(size_t l=0;l<temp.re_dim().second;l++)            out[i*b.re_dim().first+k][j*b.re_dim().second+l]=temp[k][l];    }    return out;}//由方阵上下三角部分组成的0-1阵 不包含对角线元素matrix lower_tri(const matrix&m){    if(m.re_dim().first!=m.re_dim().second)        throw runtime_error("dim error!\n");    matrix out(m.re_dim().first,m.re_dim().second,0);    for(size_t i=1;i<out.re_dim().first;i++)        for(size_t j=0;j<i;j++)        out[i][j]=1;    return out;}matrix upper_tri(const matrix&m){    return T(lower_tri(m));}//矩阵拉直(按列)c vec(const matrix&m){    c out=m.col(0);    for(size_t i=1;i<m.re_dim().second;i++)        out=connect({out,m.col(i)});    return out;}//矩阵半拉直 拉直运算忽略主对角线以上部分c vech(const matrix&m){    c out=m.col(0);    for(size_t i=1;i<m.re_dim().second;i++)    {        vector<double> temp=m.col(i).const_re_c_val();        auto a=temp.begin();        for(size_t j=0;j<i;j++)            a++;        vector<double>temp0(a,temp.end());        out=connect({out,c(temp0)});    }    return out;}//在方阵满秩情况下 给出方阵A的QR分解pair<matrix,matrix> qr_of_full_rank(const matrix&A){    if(A.re_dim().first!=A.re_dim().second)        throw runtime_error("dim error!\n");    // /*    size_t rk=rank(A);    cout<<"The rank: "<<rk<<endl;    if(rk!=A.re_dim().first)        throw runtime_error("singular!\n");    //*/    cout<<"A=QR"<<endl;    cout<<endl;    cout<<"Q of A:\n";    matrix Q=orth_normal(A);    Q.out();    cout<<"R of m0:\n";    matrix R=T(Q)%A;    R.val_to_double(lower_tri(R),0);    R.out();    // /*    cout<<"The F_norm gap between A and QR:\n";    cout<<F_norm(A-Q%R)<<endl;    //*/    return make_pair(Q,R);}//要求自变量矩阵为可逆方阵 利用qr分解 求解线性方程组//这只是一个形式上的实现 因为求解R逆的方法是用伴随矩阵法//如能改为初等行变换法应会提高效率c qr_solve(const matrix&m,const c&v){    auto qr=qr_of_full_rank(m);    return solve(qr.second)%T(qr.first)%v;}int main(){//测试程序：   matrix m2(c({5,0,10,11,0,4,4,3,10,4,3,3,11,3,3,2,10,10,10,10,11,12,13,14,15}),5,5);   m2.out();   c v1({1,1,1,1,2});   c v2=solve(m2,v1);   c v3=qr_solve(m2,v1);   v2.out();   v3.out();   cout<<"The gap of F_norm:\n";   cout<<F_norm(v2-v3)<<endl;   return 0;}