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

头文件：

[cpp] view plain copy

print?

1. class_of_R.h  
2.    
3. #ifndef CLASS_OF_R_H_  
4. #define CLASS_OF_R_H_  
5. #include<iostream>  
6. using std::cout;  
7. using std::endl;  
8. using std::cin;  
9. using std::istream;  
10. using std::ios_base;  
11. #include<vector>  
12. using std::vector;  
13. #include<initializer_list>  
14. using std::initializer_list;  
15. #include<algorithm>  
16. using std::for_each;  
17. using std::sort;  
18. using std::copy;  
19. using std::find;  
20. #include<stdexcept>  
21. using std::runtime_error;  
22. using std::out_of_range;  
23. #include<cstdlib>  
24. using std::srand;  
25. using std::rand;  
26. using std::system;  
27. #include<ctime>  
28. using std::time;  
29. #include<cmath>  
30. using std::sqrt;  
31. using std::log;  
32. using std::exp;  
33. #include<utility>  
34. using std::pair;  
35. using std::make_pair;  
36. #include<iomanip>  
37. using std::setw;  
38. #include<set>  
39. using std::multiset;  
40. using std::set;  
41. #include<unordered_map>  
42. using std::unordered_map;  
43. using std::unordered_multimap;  
44. #include<map>  
45. using std::multimap;  
46. using std::map;  
47. #include<string>  
48. using std::string;  
49. #include<fstream>  
50. using std::ifstream;  
51. #include<memory>  
52. using std::shared_ptr;  
53. #include<list>  
54. using std::list;  
55. class matrix;  
56. //类c声明  
57. class c  
58. {  
59. private:  
60.     vector<double> c_val;  
61.     size_t dim;  
62. public:  
63.     c():c_val(vector<double>()),dim(0){}  
64.     c(initializer_list<double>lst);  
65.     c(const vector<double>&c_v):c_val(c_v),dim(c_v.size()){}  
66.     c(const matrix&m);  
67.     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;}  
68.     double operator[](size_t i)const{return c_val[i];}  
69.     double& operator[](size_t i){return c_val[i];}  
70.     size_t re_dim()const{return dim;}  
71.     friend c operator+(const c&a,const c&b);  
72.     friend c operator-(const c&a);  
73.     friend c operator-(const c&a,const c&b);  
74.     friend c rep(const double&d,size_t n);  
75.     friend c operator+(const c&a,const double&d);  
76.     friend c operator-(const c&a,const double&d);  
77.     friend c operator+(const double&d,const c&a);  
78.     friend c operator-(const double&d,const c&a);  
79.     friend c operator*(const double&d,const c&a);  
80.     friend c operator*(const c&a,const double&d);  
81.     friend c operator/(const double&d,const c&a);  
82.     friend c operator/(const c&a,const double&d);  
83.     friend c operator*(const c&a,const c&b);  
84.     friend bool operator==(const c&a,const c&b){return a.c_val==b.c_val;}  
85.     friend bool operator!=(const c&a,const c&b){return a.c_val!=b.c_val;}  
86.     vector<double> const_re_c_val()const{return c_val;}  
87.     //在使用const_re_c_val()时可进行拷贝后使用  
88. };  
89. //类matrix声明  
90. class matrix  
91. {  
92. private:  
93.     vector<c> m_val;  
94.     size_t rdim;  
95.     size_t cdim;  
96. public:  
97.     matrix():m_val(vector<c>()),rdim(0),cdim(0){}  
98.     matrix(const vector<c>&m):m_val(m),rdim(m_val.size()),cdim(m_val[0].re_dim()){}  
99.     matrix(size_t m,size_t n,const double&d):m_val(vector<c>(m,vector<double>(n,d))),rdim(m),cdim(n){}  
100.     matrix(const c&v,size_t r,size_t c);  
101.     matrix(const c&v);  
102.     c operator[](size_t i)const{return m_val[i];}  
103.     c& operator[](size_t i){return m_val[i];}  
104.     void out(size_t n=3)const;  
105.     c col(size_t i)const;  
106.     void col_change(size_t i,const c&a);  
107.     friend matrix operator+(const matrix&a,const matrix&b);  
108.     friend matrix operator-(const matrix&a);  
109.     friend matrix operator-(const matrix&a,const matrix&b);  
110.     friend matrix operator+(const matrix&a,const double&d);  
111.     friend matrix operator+(const double&d,const matrix&a);  
112.     friend matrix operator-(const matrix&a,const double&d);  
113.     friend matrix operator-(const double&b,const matrix&a);  
114.     friend matrix operator*(const matrix&a,const double&d);  
115.     friend matrix operator*(const double&d,const matrix&a);  
116.     friend matrix operator/(const matrix&a,const double&d);  
117.     friend matrix operator/(const double&d,const matrix&a);  
118.     void Swap(size_t i,size_t j)  
119.     {c temp=this->operator[](i);this->operator[](i)=this->operator[](j);this->operator[](j)=temp;}  
120.     void col_Swap(size_t i,size_t j)  
121.     {  
122.         auto temp=this->col(i);  
123.         this->col_change(i,this->col(j));  
124.         this->col_change(j,temp);  
125.     }  
126.     friend bool operator==(const matrix&a,const matrix&b){return a.m_val==b.m_val;}  
127.     friend bool operator!=(const matrix&a,const matrix&b){return a.m_val!=b.m_val;}  
128.     pair<size_t,size_t> re_dim()const{return make_pair(rdim,cdim);}  
129.     friend c operator%(const matrix&m,const c&v);  
130.     friend matrix operator%(const matrix&a,const matrix&b);  
131.     friend matrix operator%(const c&v,const matrix&m);  
132. };  
133. //类factor声明  
134. class factor  
135. {  
136. private:  
137.     vector<pair<double,size_t>>f_val;  
138.     map<double,size_t>f_val_to_levels;  
139.     vector<size_t>i_val;  
140.     set<double>Levels;  
141. public:  
142.     factor(const c&v=c());  
143.     void out(size_t n=1);  
144.     friend c levels(const factor&f);  
145.     friend size_t nlevels(const factor&f);  
146.     friend factor as_integer(const factor&f);  
147.     friend factor as_factor(const c&v);  
148. };  
149.    
150. #endif // CLASS_OF_R_H_  
151.    
152. 类成员函数定义文件：  
153. class_of_R.cpp  
154.    
155. #include”class_of_R.h”  
156. //类c的有关定义  
157. c::c(initializer_list<double>lst)  
158. {  
159.     vector<double>out;  
160.     for_each(lst.begin(),lst.end(),[&](const double&d){out.push_back(d);});  
161.     c_val=out;  
162.     dim=out.size();  
163. }  
164.   
165. c::c(const matrix&m)  
166. {  
167.     vector<double>temp;  
168.   
169.         for(size_t j=0;j<m.re_dim().second;j++)  
170.                 for(size_t i=0;i<m.re_dim().first;i++)  
171.         temp.push_back(m[i][j]);  
172.     c_val=temp;  
173.     dim=temp.size();  
174. }  
175.    
176.   
177. //类matrix的有关定义  
178. void matrix::out(size_t n)const  
179. {  
180.     cout<<”The matrix is:\n”;  
181.     for(size_t i=0;i<rdim;i++)  
182.     {  
183.         for(size_t j=0;j<cdim;j++)  
184.             cout<<setw(n)<<(*this)[i][j]<<“ ”;  
185.         cout<<endl;  
186.     }  
187.     cout<<endl;  
188. }  
189. c matrix::col(size_t i)const  
190. {  
191.     vector<double>temp;  
192.     for(size_t i=0;i<rdim;i++)  
193.         temp.push_back(0);  
194.     c out(temp);  
195.     for(size_t j=0;j<rdim;j++)  
196.         out[j]=(*this)[j][i];  
197.     return out;  
198. }  
199. void matrix::col_change(size_t i,const c&a)  
200. {  
201.     for(size_t j=0;j<a.re_dim();j++)  
202.         (*this)[j][i]=a[j];  
203. }  
204. matrix::matrix(const c&v,size_t r,size_t c):rdim(r),cdim(c)  
205. {  
206.     matrix out(r,c,0);  
207.     for(size_t i=0;i<v.re_dim();i++)  
208.         out[i%r][i/r]=v[i];  
209.     (this->m_val)=out.m_val;  
210. }  
211. matrix::matrix(const c&v):rdim(v.re_dim()),cdim(1)  
212. {  
213.     matrix out(v.re_dim(),1,0);  
214.     for(size_t i=0;i<out.rdim;i++)  
215.         out[i][0]=v[i];  
216.     this->m_val=out.m_val;  
217. }  
218.    
219. //类factor有关定义  
220. factor::factor(const c&v)  
221. {  
222.     vector<double>temp=v.const_re_c_val();  
223.     for_each(temp.begin(),temp.end(),  
224.              [&](const double&d){Levels.insert(d);}  
225.              );  
226.     vector<double>s_to_v(Levels.begin(),Levels.end());  
227.     vector<pair<double,size_t>>ttemp;  
228.     for(size_t i=0;i<s_to_v.size();i++)  
229.         ttemp.push_back(make_pair(s_to_v[i],i));  
230.     for(size_t i=0;i<temp.size();i++)  
231.     {  
232.      for(size_t j=0;j<ttemp.size();j++)  
233.      {  
234.          if(temp[i]==ttemp[j].first)  
235.          {  
236.             f_val.push_back(ttemp[j]);  
237.          }  
238.       }  
239.     }  
240.   
241.     for(size_t i=0;i<ttemp.size();i++)  
242.         f_val_to_levels.insert(ttemp[i]);  
243.         auto a=f_val.begin();  
244.     for(size_t i=0;i<f_val.size();i++)  
245.     {  
246.          i_val.push_back(a->second);  
247.          a++;  
248.     }  
249. }  
250. void factor::out(size_t n)  
251. {  
252.     vector<double>temp(f_val.size(),0);  
253.     c out0(temp);  
254.     auto a=f_val.begin();  
255.     for(size_t i=0;i<out0.re_dim();i++)  
256.     {  
257.         out0[i]=a->first;  
258.         a++;  
259.     }  
260.     vector<double>temp0(Levels.size(),0);  
261.     c out1(temp0);  
262.     set<double>::iterator b=Levels.begin();  
263.     for(size_t i=0;i<out1.re_dim();i++)  
264.     {  
265.         out1[i]=*b;  
266.         b++;  
267.     }  
268.     out0.out(n);  
269.     cout<<”Levels: ”;  
270.     cout<<endl;  
271.     out1.out(n);  
272. }  
273.    
274.    
275. 其他函数即测试程序文件：  
276. use_R.cpp  
277.    
278. #include”class_of_R.h”  
279. //若干友元函数  
280. c operator+(const c&a,const c&b)  
281. {  
282.     size_t MAX=b.re_dim();  
283.     if(a.re_dim()>b.re_dim())  
284.         MAX=a.re_dim();  
285.     vector<double>temp0(MAX,0);  
286.     vector<double>temp1(MAX,0);  
287.     vector<double>temp2(MAX,0);  
288.     copy(a.c_val.begin(),a.c_val.end(),temp0.begin());  
289.     copy(b.c_val.begin(),b.c_val.end(),temp1.begin());  
290.     for(size_t i=0;i<MAX;i++)  
291.         temp2[i]=temp0[i]+temp1[i];  
292.     return c(temp2);  
293. }  
294. c operator-(const c&a)  
295. {  
296.     c out(a);  
297.     for(size_t i=0;i<out.re_dim();i++)  
298.         out[i]=-out[i];  
299.     return out;  
300. }  
301. c operator-(const c&a,const c&b)  
302. {  
303.     return (a+(-b));  
304. }  
305. c rep(const double&d,size_t n)  
306. {  
307.     return c(vector<double>(n,d));  
308. }  
309. c operator+(const c&a,const double&d)  
310. {  
311.     return a+rep(d,a.re_dim());  
312. }  
313. c operator-(const c&a,const double&d)  
314. {  
315.     return (a+(-d));  
316. }  
317. c operator+(const double&d,const c&a)  
318. {  
319.     return a+d;  
320. }  
321. c operator-(const double&d,const c&a)  
322. {  
323.     return (-a)+d;  
324. }  
325. c operator*(const double&d,const c&a)  
326. {  
327.     c out(a);  
328.     for(size_t i=0;i<a.re_dim();i++)  
329.         out[i]*=d;  
330.     return out;  
331. }  
332. c operator*(const c&a,const double&d)  
333. {  
334.     return d*a;  
335. }  
336. c operator/(const double&d,const c&a)  
337. {  
338.     c out(a);  
339.     for(size_t i=0;i<a.re_dim();i++)  
340.         out[i]=d/a[i];  
341.     return out;  
342. }  
343. c operator/(const c&a,const double&d)  
344. {  
345.     c out(a);  
346.     for(size_t i=0;i<a.re_dim();i++)  
347.         out[i]=a[i]/d;  
348.     return out;  
349. }  
350. c operator*(const c&a,const c&b)  
351. {  
352.     if(a.re_dim()!=b.re_dim())  
353.         throw runtime_error(“dim error!\n”);  
354.     c out(a);  
355.     for(size_t i=0;i<a.re_dim();i++)  
356.         out[i]*=b[i];  
357.     return out;  
358. }  
359. matrix operator+(const matrix&a,const matrix&b)  
360. {  
361.     size_t r_MAX=a.rdim,c_MAX=a.cdim;  
362.     if(a.rdim<b.rdim)  
363.         r_MAX=b.rdim;  
364.     if(a.cdim<b.cdim)  
365.         c_MAX=b.cdim;  
366.     matrix tempa(r_MAX,c_MAX,0);  
367.     for(size_t i=0;i<a.rdim;i++)  
368.         for(size_t j=0;j<a.cdim;j++)  
369.         tempa[i][j]=a[i][j];  
370.     matrix tempb(r_MAX,c_MAX,0);  
371.     for(size_t i=0;i<b.rdim;i++)  
372.         for(size_t j=0;j<b.cdim;j++)  
373.         tempb[i][j]=b[i][j];  
374.     matrix tempout(r_MAX,c_MAX,0);  
375.     for(size_t i=0;i<tempout.rdim;i++)  
376.         for(size_t j=0;j<tempout.cdim;j++)  
377.         tempout[i][j]=tempa[i][j]+tempb[i][j];  
378.     return tempout;  
379. }  
380. matrix operator-(const matrix&a)  
381. {  
382.     matrix out(a);  
383.     for(size_t i=0;i<out.rdim;i++)  
384.         for(size_t j=0;j<out.cdim;j++)  
385.         out[i][j]=-out[i][j];  
386.     return out;  
387. }  
388. matrix operator-(const matrix&a,const matrix&b)  
389. {  
390.     return ((-b)+a);  
391. }  
392. matrix operator+(const matrix&a,const double&d)  
393. {  
394.     return a+matrix(a.rdim,a.cdim,d);  
395. }  
396. matrix operator+(const double&d,const matrix&a)  
397. {  
398.     return a+d;  
399. }  
400. matrix operator-(const matrix&a,const double&d)  
401. {  
402.     return a+(-d);  
403. }  
404. matrix operator-(const double&b,const matrix&a)  
405. {  
406.     return (-a)+b;  
407. }  
408. matrix operator*(const matrix&a,const double&d)  
409. {  
410.     matrix out(a);  
411.     for(size_t i=0;i<out.rdim;i++)  
412.         for(size_t j=0;j<out.cdim;j++)  
413.         out[i][j]*=d;  
414.     return out;  
415. }  
416. matrix operator*(const double&d,const matrix&a)  
417. {  
418.     return a*d;  
419. }  
420. matrix operator/(const matrix&a,const double&d)  
421. {  
422.     return a*(1/d);  
423. }  
424. matrix operator/(const double&d,const matrix&a)  
425. {  
426.     matrix out(a.rdim,a.cdim,0);  
427.     for(size_t i=0;i<out.rdim;i++)  
428.         for(size_t j=0;j<out.cdim;j++)  
429.         out[i][j]=d/a[i][j];  
430.     return out;  
431. }  
432. c levels(const factor&f)  
433. {  
434.     vector<double> temp(f.Levels.begin(),f.Levels.end());  
435.     c done(temp);  
436.     done.out();  
437.     return done;  
438. }  
439. size_t nlevels(const factor&f)  
440. {  
441.     size_t out=f.Levels.size();  
442.     cout<<out<<endl;  
443.     return out;  
444. }  
445. //返回输入因子相应在存储中的下标构成的因子  
446. factor as_integer(const factor&f)  
447. {  
448.     vector<double> temp(f.i_val.begin(),f.i_val.end());  
449.     factor temp0=factor(c(temp));  
450.     for_each(temp0.f_val.begin(),temp0.f_val.end(),  
451.              [](const pair<double,size_t>&p)  
452.              {  
453.                  cout<<p.first<<” ”;  
454.              }  
455.              );  
456.              cout<<endl;  
457.     return temp0;  
458. }  
459. factor as_factor(const c&v)  
460. {  
461.     factor temp0=factor(v);  
462.     temp0.out();  
463.     return temp0;  
464. }  
465. //基本运算函数  
466. matrix diag(const c&v)  
467. {  
468.     matrix out(v.re_dim(),v.re_dim(),0);  
469.     for(size_t i=0;i<v.re_dim();i++)  
470.         out[i][i]=v[i];  
471.     return out;  
472. }  
473. matrix diag(const matrix&m)  
474. {  
475.     if(m.re_dim().first!=m.re_dim().second)  
476.         throw runtime_error(“dim error!\n”);  
477.     matrix out(m.re_dim().first,m.re_dim().second,0);  
478.     for(size_t i=0;i<out.re_dim().first;i++)  
479.         out[i][i]=m[i][i];  
480.     return out;  
481. }  
482. double inter_product(const c&a,const c&b)  
483. {  
484.     if(a.re_dim()!=b.re_dim())  
485.         throw runtime_error(“dim error!\n”);  
486.     double out=0;  
487.     for(size_t i=0;i<a.re_dim();i++)  
488.     out+=(a[i]*b[i]);  
489.     return out;  
490. }  
491. c operator%(const matrix&m,const c&v)  
492. {  
493.     if(m.re_dim().second!=v.re_dim())  
494.         throw runtime_error(“dim error!\n”);  
495.     c out(rep(0,m.re_dim().first));  
496.     for(size_t i=0;i<out.re_dim();i++)  
497.         out[i]=inter_product(m[i],v);  
498.     return out;  
499. }  
500. //同下  
501. double prod(const c&v)  
502. {  
503.     double out=1;  
504.     for(size_t i=0;i<v.re_dim();i++)  
505.         out*=v[i];  
506.     return out;  
507. }  
508. //对于matrix的求和由c复制构造函数的重载函数定义：c(const matrix&m)  
509. double sum(const c&v)  
510. {  
511.     double out=v[0];  
512.     for(size_t i=1;i<v.re_dim();i++)  
513.         out+=v[i];  
514.     return out;  
515. }  
516. matrix operator%(const c&v,const matrix&m)  
517. {  
518.     if(v.re_dim()!=m.re_dim().first)  
519.         throw runtime_error(“dim error!\n”);  
520.     matrix out(1,m.re_dim().second,0);  
521.     for(size_t i=0;i<out.re_dim().second;i++)  
522.         out[0][i]=inter_product(v,m.col(i));  
523.     return out;  
524. }  
525. matrix outer_product(const c&a,const c&b)  
526. {  
527.     matrix out(a.re_dim(),b.re_dim(),0);  
528.     for(size_t i=0;i<out.re_dim().first;i++)  
529.         for(size_t j=0;j<out.re_dim().second;j++)  
530.         out[i][j]=a[i]*b[j];  
531.     return out;  
532. }  
533. size_t length(const c&v)  
534. {  
535.     return v.re_dim();  
536. }  
537. c sort(const c&v)  
538. {  
539.     vector<double>out;  
540.     for(size_t i=0;i<v.re_dim();i++)  
541.         out.push_back(v[i]);  
542.     sort(out.begin(),out.end());  
543.     return c(out);  
544. }  
545. double mean(const c&v)  
546. {  
547.     return sum(v)/length(v);  
548. }  
549. //其转换版本返回矩阵转为向量的方差  
550. double var(const c&v)  
551. {  
552.     return mean((v-mean(v))*(v-mean(v)));  
553. }  
554. c connect(initializer_list<c> lst)  
555. {  
556.     vector<double>temp;  
557.     for_each(lst.begin(),lst.end(),[&](const c&v){for(auto a:v.const_re_c_val())temp.push_back(a);});  
558.     return c(temp);  
559. }  
560. //标准差  
561. template<typename any>  
562. double STD(const any&a)  
563. {  
564.     return sqrt(var(a)*(length(a))/(length(a)-1));  
565. }  
566. matrix rbind(initializer_list<c>lst)  
567. {  
568.     size_t r=0;  
569.     auto a=lst.begin();  
570.     while(a!=lst.end())  
571.     {  
572.         a++;  
573.         r++;  
574.     }  
575.     size_t co=0;  
576.     for_each(lst.begin(),lst.end(),[&](const c&v){if(v.re_dim()>co)co=v.re_dim();});  
577.     matrix out(r,co,0);  
578.     size_t i=0;  
579.     a=lst.begin();  
580.     while(a!=lst.end())  
581.     {  
582.         vector<double>temp(co,0);  
583.         for(size_t i=0;i<a->const_re_c_val().size();i++)  
584.             temp[i]=a->operator[](i);  
585.         out[i]=c(temp);  
586.         i++;  
587.         a++;  
588.     }  
589.   
590.     return out;  
591. }  
592. size_t dim(const c&v)  
593. {  
594.     return v.re_dim();  
595. }  
596. pair<size_t,size_t> dim(const matrix&m)  
597. {  
598.     return m.re_dim();  
599. }  
600.   
601. double min(const c&v)  
602. {  
603.     double temp=v[0];  
604.     for(size_t i=1;i<length(v);i++)  
605.         if(v[i]<temp)  
606.         temp=v[i];  
607.     return temp;  
608. }  
609. double max(const c&v)  
610. {  
611.     double temp=v[0];  
612.     for(size_t i=1;i<length(v);i++)  
613.         if(v[i]>temp)  
614.         temp=v[i];  
615.     return temp;  
616. }  
617. c pmin(initializer_list<c>lst)  
618. {  
619.     c temp=rep(0,lst.size());  
620.     auto a=lst.begin();  
621.     for(size_t i=0;i<temp.re_dim();i++)  
622.     {  
623.         temp[i]=min(*a);  
624.         a++;  
625.     }  
626.     return temp;  
627. }  
628. c pmax(initializer_list<c>lst)  
629. {  
630.     c temp=rep(0,lst.size());  
631.     auto a=lst.begin();  
632.     for(size_t i=0;i<temp.re_dim();i++)  
633.     {  
634.         temp[i]=max(*a);  
635.         a++;  
636.     }  
637.     return temp;  
638. }  
639. template<typename any>  
640. pair<double,double> range(const any&v)  
641. {  
642.     cout<<”MIN: ”<<min(v)<<“ MAX: ”<<max(v)<<endl;  
643.     return make_pair(min(v),max(v));  
644. }  
645. double median(const c&v)  
646. {  
647.     return sort(v)[length(v)/2];  
648. }  
649. double var(const c&a,const c&b)  
650. {  
651.     if(length(a)!=length(b))  
652.         throw runtime_error(“dim error!\n”);  
653.     c t1(a-mean(a)),t2(b-mean(b));  
654.     return inter_product(t1,t2)/(length(a)-1);  
655. }  
656. double cov(const c&a,const c&b)  
657. {  
658.     return var(a,b);  
659. }  
660. double cor(const c&a,const c&b)  
661. {  
662.     return var(a,b)/STD(a)/STD(b);  
663. }  
664. matrix var(const matrix&a,const matrix&b)  
665. {  
666.     if(a.re_dim().first!=b.re_dim().first)  
667.         throw runtime_error(“dim error!\n”);  
668.     matrix temp(a.re_dim().second,b.re_dim().second,0);  
669.     for(size_t i=0;i<temp.re_dim().first;i++)  
670.         for(size_t j=0;j<temp.re_dim().second;j++)  
671.         temp[i][j]=var(a.col(i),b.col(j));  
672.     return temp;  
673. }  
674. matrix cov(const matrix&a,const matrix&b)  
675. {  
676.     return var(a,b);  
677. }  
678. matrix cor(const matrix&a,const matrix&b)  
679. {  
680.     matrix temp=var(a,b);  
681.     for(size_t i=0;i<temp.re_dim().first;i++)  
682.         for(size_t j=0;j<temp.re_dim().second;j++)  
683.         temp[i][j]/=((STD(a.col(i))*STD(b.col(j))));  
684.     return temp;  
685. }  
686. matrix cov(const matrix&a)  
687. {  
688.     return var(a,a);  
689. }  
690. matrix cor(const matrix&a)  
691. {  
692.     return cor(a,a);  
693. }  
694. //由于关联容器的删除方法利用误差排序的方法  
695. //利用顺序容器迭代器的删除方法有困难  
696. c rank(const c&v)  
697. {  
698.     vector<double>temp00(v.const_re_c_val());  
699.     for(size_t i=0;i<temp00.size();i++)  
700.         temp00[i]+=(0.00000001*i);  
701.     vector<double>temp0(temp00);  
702.     multiset<double>temp2(temp0.begin(),temp0.end());  
703.     vector<double>temp1(v.re_dim());  
704.     size_t i=0,j=0;  
705.     double pre_val=min(c(vector<double>(temp2.begin(),temp2.end())));  
706.     vector<double>::iterator s0;  
707.     while(i<v.re_dim())  
708.     {  
709.     double val=min(c(vector<double>(temp2.begin(),temp2.end())));  
710.     s0=find(temp0.begin(),temp0.end(),val);  
711.     if(sqrt((pre_val-val)*(pre_val-val))>0.000001)  
712.         j++;  
713.     temp1[s0-temp0.begin()]=j;  
714.     temp2.erase(val);  
715.     i++;  
716.     pre_val=val;  
717.     }  
718.     return  c(temp1);  
719. }  
720. c rev(const c&v)  
721. {  
722.     vector<double> temp;  
723.     auto a=v.const_re_c_val().rbegin();  
724.     while(a!=v.const_re_c_val().rend())  
725.     {  
726.         temp.push_back(*a);  
727.         a++;  
728.     }  
729.     return c(temp);  
730. }  
731. //向量先排序再输出排序前的下标值  
732. c order(const c&v)  
733. {  
734.     unordered_multimap<double,double>m0;  
735.     for(size_t i=0;i<length(v);i++)  
736.         m0.insert(make_pair(v[i],i));  
737.     multimap<double,double>m1(m0.begin(),m0.end());  
738.     c temp(rep(0,length(v)));  
739.     auto a=m1.begin();  
740.     size_t i=0;  
741.     while(a!=m1.end())  
742.     {  
743.         temp[i]=a->second;  
744.         i++;  
745.         a++;  
746.     }  
747.     return temp;  
748. }  
749. c cumsum(const c&v)  
750. {  
751.     c temp(v);  
752.     for(size_t i=1;i<dim(v);i++)  
753.         temp[i]+=temp[i-1];  
754.     return temp;  
755. }  
756. c cumprod(const c&v)  
757. {  
758.     c temp(v);  
759.     for(size_t i=1;i<length(v);i++)  
760.         temp[i]*=temp[i-1];  
761.     return temp;  
762. }  
763. matrix cumsum(const matrix&m)  
764. {  
765.     c temp(m);  
766.     for(size_t i=1;i<length(m);i++)  
767.         temp[i]+=temp[i-1];  
768.     return matrix(temp,m.re_dim().first,m.re_dim().second);  
769. }  
770. matrix cumprod(const matrix&m)  
771. {  
772.     c temp(m);  
773.     for(size_t i=1;i<length(m);i++)  
774.         temp[i]*=temp[i-1];  
775.     return matrix(temp,m.re_dim().first,m.re_dim().second);  
776. }  
777.    
778. //矩阵运算函数  
779. //初等行变换化上三角求解行列式  
780. double det(const matrix&m)  
781. {  
782. auto f=[](const matrix&mm)->pair<matrix,double>  
783. {  
784.     if(mm.re_dim().first!=mm.re_dim().second)  
785.        throw runtime_error(“dim error!\n”);  
786.     size_t MIN=mm.re_dim().first;  
787.   
788.     double mult=1;  
789.     size_t col=0;  
790.     size_t bar=0;  
791.     size_t b_c=0;  
792.     matrix t(mm);  
793.   
794.     while(b_c<MIN)  
795.     {  
796.     for(size_t i=bar+1;i<MIN;)  
797.     {  
798.         if(t[i][col]==0)  
799.             i++;  
800.         else  
801.         {  
802.            t.Swap(i,bar);  
803.            mult*=(-1);  
804.            break;  
805.         }  
806.     }  
807.     if(bar!=MIN+1)  
808.     {  
809.     if(t[bar][col]!=0)  
810.     {  
811.        for(size_t i=bar+1;i<MIN;i++)  
812.      {  
813.         t[bar][i]=t[bar][i]/t[bar][col];  
814.         if(i==bar+1)  
815.         mult*=t[bar][col];  
816.      }  
817.         if(col!=MIN-1)  
818.         t[bar][col]=1;  
819.     }  
820.   
821.      for(size_t i=bar+1;i<MIN;i++)  
822.      {  
823.          t[i]=t[i]-t[i][col]*t[bar];  
824.      }  
825.     }  
826.       col++;  
827.       bar++;  
828.       b_c=bar;  
829.       if(col>bar);  
830.       b_c=col;  
831.     }  
832. return make_pair(t,mult);  
833. };  
834.     pair<matrix,double> tt=f(m);  
835.     return tt.second*prod(diag(tt.first)%rep(1,m.re_dim().first));  
836. }  
837. //Cramer法则解线性方程组：要求自变量矩阵为方阵且相应行列式不为0  
838. c solve(const matrix&m,const c&v)  
839. {  
840.     c temp(rep(0,v.re_dim()));  
841.     double low=det(m);  
842.     if(m.re_dim().first!=m.re_dim().second)  
843.       throw runtime_error(“dim error!\n”);  
844.     if(low==0)  
845.       throw runtime_error(“singular!\n”);  
846.     for(size_t i=0;i<v.re_dim();i++)  
847.     {  
848.         matrix temp0(m);  
849.         temp0.col_change(i,v);  
850.         temp[i]=det(temp0)/low;  
851.     }  
852.     return temp;  
853. }  
854. double tr(const matrix&m)  
855. {  
856.     if(m.re_dim().first!=m.re_dim().second)  
857.         throw runtime_error(“dim error!\n”);  
858.     return inter_product(diag(m)%rep(1,m.re_dim().first),rep(1,m.re_dim().first));  
859. }  
860. matrix T(const matrix&m)  
861. {  
862.     matrix out(m.re_dim().second,m.re_dim().first,0);  
863.     for(size_t i=0;i<out.re_dim().first;i++)  
864.         for(size_t j=0;j<out.re_dim().second;j++)  
865.         out[i][j]=m[j][i];  
866.     return out;  
867. }  
868. matrix cbind(initializer_list<c>lst)  
869. {  
870.     return T(rbind(lst));  
871. }  
872. matrix operator%(const matrix&a,const matrix&b)  
873. {  
874.     matrix out(a.re_dim().first,b.re_dim().second,0);  
875.     for(size_t i=0;i<out.re_dim().first;i++)  
876.         for(size_t j=0;j<out.re_dim().second;j++)  
877.         out[i][j]=inter_product(a[i],b.col(j));  
878.     return out;  
879. }  
880. double F_norm(const matrix&m)  
881. {  
882.     return sqrt(tr(T(m)%m));  
883. }  
884.    
885. //a,b为区间 d为步长 e为精度 默认为求整数精度：以1为精度特征值 如所求特征值太少可通过减少d改善  
886. //只能求实值  
887. //可以用F_norm作为a、b的界（-F_norm(m),F_norm(m)）可确保值在范围中  
888. //可用于求解但效率低  
889. c eigen(const matrix&m,const double&a,const double &b,const double &d=0.0001,const double &e=1)  
890. {  
891.     if(m.re_dim().first!=m.re_dim().second)  
892.         throw runtime_error(“dim error!\n”);  
893.     vector<double>temp;  
894.     double pre_beg=100000000;  
895.     for(double beg=a;beg<b;beg+=d)  
896.     {  
897.         matrix mm0(diag(rep(beg,m.re_dim().first))-m);  
898.         double s0=det(mm0);  
899.         if(sqrt(s0*s0)<e)  
900.         {  
901.             if(sqrt((pre_beg-beg)*(pre_beg-beg))>1)  
902.             {  
903.              temp.push_back(beg);  
904.              pre_beg=beg;  
905.              cout<<”Done!\n”;  
906.              cout<<beg<<endl;//eigenvalue  
907.             }  
908.         }  
909.     }  
910.     return c(temp);  
911. }  
912. //重载求eign的函数 d:步长 可保证绝对求出所有特征值（利用循环增大精度）  
913. c eigen(const matrix&m,const double&d=0.1)  
914. {  
915.     double dd=d;  
916.     double ra=F_norm(m);  
917.     double trr=tr(m);  
918.     c temp;  
919.     double r=0;  
920.     do  
921.     {  
922.         dd*=0.1;//原来设为0.1  
923.         temp=eigen(m,-ra,ra,dd);  
924.         if(temp.re_dim()==m.re_dim().first)  
925.             break;  
926.         r=sum(temp)-trr;  
927.     }  
928.     while(sqrt((r*r)>1));  
929.         //与迹的差距（精度）设定为1 绝对值小于此值的特征值被忽略  
930.         //此种方法的一个缺陷为异号特征值的抵消 对于同号特征值有把握  
931.         //只适用于所有特征值为实值的情况  
932.         //如有复特征值，只输出实数的并循环不能结束  
933.     return temp;  
934. }  
935. //矩阵的余子阵  
936. matrix M(const matrix&m,size_t k,size_t l)  
937. {  
938.     matrix temp(m);  
939.     set<size_t> s1;  
940.     for(size_t i=0;i<m.re_dim().first;i++)  
941.         s1.insert(i);  
942.     set<size_t> s2;  
943.     for(size_t i=0;i<m.re_dim().second;i++)  
944.         s2.insert(i);  
945.     s1.erase(k);  
946.     s2.erase(l);  
947.     vector<double>temp00;  
948.     for(auto j:s2)  
949.         for(auto i:s1)  
950.         temp00.push_back(temp[i][j]);  
951.     return matrix (c(temp00),m.re_dim().first-1,m.re_dim().second-1);  
952. }  
953. //方阵的伴随矩阵  
954. matrix A_accompany(const matrix&m)  
955. {  
956.     matrix out(m.re_dim().first,m.re_dim().second,0);  
957.     for(size_t i=0;i<out.re_dim().first;i++)  
958.         for(size_t j=0;j<out.re_dim().second;j++)  
959.         out[i][j]=det(M(m,i,j))*(((i+j)%2)?-1:1);  
960.         //代数余子式  
961.     return T(out);  
962. }  
963. //求方阵的逆矩阵（利用伴随矩阵方法）  
964. matrix solve(const matrix&m)  
965. {  
966.     double de=det(m);  
967.     if(sqrt(de*de)<0.000001)  
968.      throw runtime_error(“singular!\n”);  
969.     return A_accompany(m)/de;  
970. }  
971. //利用逆矩阵法求解方程组 要求自变量阵为方阵 且可逆  
972. c solve_eq(const matrix&m,const c&v)  
973. {  
974.     c out;  
975.     double de=det(m);  
976.     if(sqrt(de*de)<0.000001)  
977.       throw runtime_error(“singular!\n”);  
978.     return solve(m)%v;  
979. }  
980. //化列向量组构成的矩阵为正交化后列向量组构成的矩阵（Schmidt 正交化）  
981. matrix orth_normal(const matrix &m)  
982. {  
983.     matrix out(rep(0,m.re_dim().first*m.re_dim().second),m.re_dim().first,m.re_dim().second);  
984.     vector<c>dec;  
985.     dec.push_back(m.col(0)/sqrt(inter_product(m.col(0),m.col(0))));  
986.   
987.     for(size_t i=1;i<m.re_dim().second;i++)  
988.     {  
989.         c temp=m.col(i);  
990.         for(size_t j=1;j<=i;j++)  
991.         temp=temp-inter_product(dec[j-1],m.col(i))*dec[j-1];  
992.         temp=temp/sqrt(inter_product(temp,temp));  
993.         dec.push_back(temp);  
994.     }  
995.     for(size_t i=0;i<m.re_dim().second;i++)  
996.         out.col_change(i,dec[i]);  
997.     return out;  
998. }  
999. //生成相应位数的由0/1字符构成的字符串  
1000. vector<string> generator0_1(size_t n)  
1001. {  
1002.     list<string>out;  
1003.     list<string>require;  
1004.     out.push_back(”0”);  
1005.     out.push_back(”1”);  
1006.     size_t i=0;  
1007.     while(i<n-1)  
1008.     {  
1009.         if(i>0)  
1010.          out=require;  
1011.          size_t s0=out.size();  
1012.     list<string>::iterator pre_begin=out.begin();  
1013.     auto pre_end=out.end();  
1014.     size_t k=0;  
1015.     while(pre_begin!=pre_end)  
1016.     {  
1017.         pre_begin++;  
1018.         k++;  
1019.     }  
1020.     pre_begin=out.begin();  
1021.    
1022.         for(size_t i=0;i<k;i++)  
1023.       {  
1024.           string temp0=*pre_begin+”0”;  
1025.           string temp1=*pre_begin+”1”;  
1026.          out.push_back(temp0);  
1027.          out.push_back(temp1);  
1028.          pre_begin++;  
1029.       }  
1030.       require.clear();  
1031.       list<string>::iterator a=out.begin();  
1032.       for(size_t i=0;i<s0;i++)  
1033.         a++;  
1034.       while(a!=out.end())  
1035.       {  
1036.           require.push_back(*a);  
1037.           a++;  
1038.       }  
1039. i++;  
1040.     }  
1041.    return vector<string>(require.begin(),require.end());  
1042. }  
1043. //利用矩阵向量方法生成上述数组  
1044. matrix matrix_generator0_1(size_t n)  
1045. {  
1046.     matrix temp0(2,1,0);  
1047.     matrix temp1;  
1048.     temp0[0][0]=0;  
1049.     temp0[1][0]=1;  
1050.     size_t i=0;  
1051.   
1052.     while(i<n-1)  
1053.     {  
1054.         if(i>0)  
1055.          temp0=temp1;  
1056. //         size_t s0=temp0.re_dim().second;  
1057.     size_t pre_begin=0;  
1058.     auto pre_end=temp0.re_dim().first;  
1059.     temp1=matrix(temp0.re_dim().first*2,temp0.re_dim().second+1,0);  
1060.     size_t temp1_index=0;  
1061.         while(pre_begin!=pre_end)  
1062.       {  
1063.          vector<double> tt0=temp0[pre_begin].const_re_c_val();  
1064.          vector<double> tt1=tt0;  
1065.          tt0.push_back(0);  
1066.          tt1.push_back(1);  
1067.          temp1[temp1_index]=c(tt0);  
1068.          temp1_index++;  
1069.          temp1[temp1_index]=c(tt1);  
1070.          temp1_index++;  
1071.          pre_begin++;  
1072.       }  
1073.       i++;  
1074.     }  
1075.   
1076.    return temp1;  
1077. }  
1078.   
1079. //将固定位数的下标值与0/1数值构成的map以vector的形式返回  
1080. //每个map都是一种组合  
1081. vector<map<size_t,size_t>> index_map_generator(size_t n)  
1082. {  
1083.     vector<string>s0(generator0_1(n));  
1084.     vector<map<size_t,size_t>>outer;  
1085.     for(size_t i=0;i<s0.size();i++)  
1086.     {  
1087.         map<size_t,size_t>tp;  
1088.         for(size_t j=0;j<s0[i].size();j++)  
1089.        {  
1090.         size_t temp=1;  
1091.         if(s0[i][j]==‘0’)  
1092.         temp=0;  
1093.         tp.insert(make_pair(j,temp));  
1094.        }  
1095.         outer.push_back(tp);  
1096.     }  
1097.     return outer;  
1098. }  
1099. //上述产生方式的矩阵版本  
1100. vector<map<size_t,size_t>> index_map_matrix_generator(size_t n)  
1101. {  
1102.     matrix temp=matrix_generator0_1(n);  
1103.     vector<map<size_t,size_t>> out(temp.re_dim().first);  
1104.     for(size_t i=0;i<out.size();i++)  
1105.     {  
1106.         map<size_t,size_t>m;  
1107.         for(size_t j=0;j<temp.re_dim().second;j++)  
1108.             m.insert(make_pair(j,temp[i][j]));  
1109.         out[i]=m;  
1110.     }  
1111.     return out;  
1112. }  
1113. //返回上述index_map_generator(m)中指向0_1和为n的子集  
1114. vector<map<size_t,size_t>> require_dim(size_t m,size_t n)  
1115. {  
1116.     vector<map<size_t,size_t>>v0=index_map_generator(m);  
1117.     vector<map<size_t,size_t>>out;  
1118.     auto a=v0.begin();  
1119.     while(a!=v0.end())  
1120.     {  
1121.         size_t temp=0;  
1122.         for_each(a->begin(),a->end(),[&](const pair<size_t,size_t>&p){  
1123.                  if(p.second==1)  
1124.                     temp++;  
1125.                  });  
1126.         if(temp==n)  
1127.             out.push_back(*a);  
1128.         a++;  
1129.     }  
1130.     return out;  
1131. }  
1132. //matrix 所有n阶子式构成的向量  
1133. vector<double> seq_of_subdet(const matrix&m,size_t n)  
1134. {  
1135.     vector<double>out;  
1136.     vector<map<size_t,size_t>>r_map=require_dim(m.re_dim().first,n);  
1137.     vector<map<size_t,size_t>>c_map=require_dim(m.re_dim().second,n);  
1138.     for(size_t i=0;i<r_map.size();i++)  
1139.         for(size_t j=0;j<c_map.size();j++)  
1140.     {  
1141.         set<size_t>r_require;  
1142.         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);});  
1143.         set<size_t>c_require;  
1144.         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);});  
1145.         matrix r_matrix(n,m.re_dim().second,0);  
1146.         auto a=r_require.begin();  
1147.         for(size_t i=0;i<r_matrix.re_dim().first;i++)  
1148.         {  
1149.             r_matrix[i]=m[*a];  
1150.             a++;  
1151.         }  
1152.         matrix c_matrix(r_matrix.re_dim().first,n,0);  
1153.         auto b=c_require.begin();  
1154.         for(size_t i=0;i<n;i++)  
1155.         {  
1156.             c_matrix.col_change(i,r_matrix.col(*b));  
1157.             b++;  
1158.         }  
1159.         out.push_back(det(c_matrix));  
1160.     }  
1161.     return out;  
1162. }  
1163. //利用子式的方法求出矩阵的秩  
1164. size_t rank(const matrix&m)  
1165. {  
1166.     size_t MAX=m.re_dim().first;  
1167.     if(m.re_dim().second<m.re_dim().first)  
1168.         MAX=m.re_dim().second;  
1169.     size_t out=MAX;  
1170.     size_t pre_num_of_non0=0;  
1171.     for(size_t i=MAX;i>0;)  
1172.     {  
1173.         vector<double>temp=seq_of_subdet(m,i);  
1174.         auto a=temp.begin();  
1175.         while(a!=temp.end())  
1176.         {  
1177.             *a=sqrt((*a)*(*a));  
1178.             a++;  
1179.         }  
1180.         size_t num_of_non0=0;  
1181.         for_each(temp.begin(),temp.end(),[&](const double&d){if(d>0.00000001)num_of_non0++;});  
1182.         if(pre_num_of_non0==0)  
1183.             i–;  
1184.         if(num_of_non0!=0)  
1185.             return i+1;  
1186.         pre_num_of_non0=num_of_non0;  
1187.     }  
1188.     return out;  
1189. }  
1190. //回归系数的求解 要求增广后的数据阵（含常系数）可逆  
1191. c lm(const c&y,const matrix&m)  
1192. {  
1193.     matrix X=matrix(m.re_dim().first,m.re_dim().second+1,0);  
1194.     X.col_change(0,rep(1,X.re_dim().first));  
1195.     for(size_t i=0;i<m.re_dim().second;i++)  
1196.         X.col_change(i+1,m.col(i));  
1197.     if(det(T(X)%X)<0.00000001)  
1198.         throw runtime_error(“singular!\n”);  
1199.     c out(solve(T(X)%X)%T(X)%y);  
1200.     cout<<”The intercept is: ”<<out[0]<<endl;  
1201.     cout<<”Others are:\n”;  
1202.     for(size_t i=1;i<out.re_dim();i++)  
1203.         cout<<”x”<<i<<“ is ”<<out[i]<<endl;  
1204.     cout<<endl;  
1205.     return out;  
1206. }  
1207. //抽样函数 参数为向量v 抽取的个数n 是否有放回replace 抽取的概率prob  
1208. //使用此方法前必须在函数调用外设定随机数种子  
1209. //要求没有0概率点  
1210. c sample(const c&v,size_t n,bool replace=1,const c&prob=c())  
1211. {  
1212.     if(!replace&&v.re_dim()<n)  
1213.         throw out_of_range(“n is too large\n”);  
1214.         //抽样概率  
1215.     vector<double>probb;  
1216.     if(prob==c())  
1217.         probb=rep(1/(double)v.re_dim(),v.re_dim()).const_re_c_val();  
1218.     else  
1219.     {  
1220.         probb=prob.const_re_c_val();  
1221.         double total=0;  
1222.         for(auto a:probb)  
1223.             total+=a;  
1224.         for(auto &a:probb)  
1225.             a/=total;  
1226.     }  
1227.     ifstream fin(”h://runif0.txt”);  
1228.     vector<double>runif;  
1229.     double val;  
1230.     while(fin>>val)  
1231.     {  
1232.         runif.push_back(val);  
1233.     }  
1234.     fin.clear();  
1235.   
1236.     //c(probb).out();  
1237.     c cprobb=cumsum(c(probb));  
1238.     //cprobb.out();  
1239.     c require(rep(0,n));  
1240.     if(replace)  
1241.     {  
1242.           size_t i=0;  
1243.           while(i<n)  
1244.           {  
1245.           double uniform=runif[rand()%30000];  
1246.           //cout<<uniform<<endl;  
1247.           double conclusion=v[v.re_dim()-1];  
1248.           for(size_t i=0;cprobb[i]<1;i++)  
1249.           {  
1250.               if(uniform<cprobb[i])  
1251.               {  
1252.                 conclusion=v[i];  
1253.                 break;  
1254.               }  
1255.           }  
1256.           require[i]=conclusion;  
1257.           i++;  
1258.           }  
1259.     }  
1260.     else  
1261.         {  
1262.             size_t i=0;  
1263.             map<double,double>mdd;  
1264.             for(size_t i=0;i<v.re_dim();i++)  
1265.                 mdd.insert(make_pair(v[i],c(probb)[i]));  
1266.             while(i<n)  
1267.             {  
1268.             vector<double> temp_v;  
1269.             vector<double> temp_prob;  
1270.             for_each(mdd.begin(),mdd.end(),[&](const pair<double,double>&p){  
1271.                     temp_v.push_back(p.first);  
1272.                     temp_prob.push_back(p.second);  
1273.                      });  
1274.             double d=sample(c(temp_v),1,1,c(temp_prob)).const_re_c_val()[0];  
1275.             require[i]=d;  
1276.             mdd.erase(d);  
1277.             i++;  
1278.             }  
1279.         }  
1280.   
1281. return require;  
1282. }  
1283. //产生几何分布随机数 n 产生个数 p成功概率  
1284. c rgeom(const size_t& n,const double& p)  
1285. {  
1286.     vector<double>out;  
1287.     size_t i=0;  
1288.     srand(0);  
1289.     ifstream fin(”h://runif0.txt”);  
1290.     vector<double>runif;  
1291.     double val;  
1292.     while(fin>>val)  
1293.     {  
1294.         runif.push_back(val);  
1295.     }  
1296.     fin.clear();  
1297.     while(i<n)  
1298.     {  
1299.         double temp=runif[rand()%30000];  
1300.         out.push_back(size_t(log(temp)/log(1-p))+1);  
1301.         i++;  
1302.     }  
1303.     return c(out);  
1304. }  
1305. //possion分布随机数 n 生成个数 d 分布参数  
1306. c rpois(const size_t&n,const double&d)  
1307. {  
1308.     size_t i=0;  
1309.     vector<double>conclusion;  
1310.     srand(time(0));  
1311.     ifstream fin(”h://runif0.txt”);  
1312.     vector<double>runif;//储存30000个均匀随机数的vector  
1313.     double val;  
1314.     while(fin>>val)  
1315.     {  
1316.         runif.push_back(val);  
1317.     }  
1318.     fin.clear();  
1319.   
1320.     while(i<n)  
1321.     {  
1322.         double p=exp(-1*d);  
1323.         size_t out=0;  
1324.         double F=p;  
1325.         double temp=runif[rand()%30000];  
1326.   
1327.          while(temp>=F)  
1328.          {  
1329.             p=d/(double)(out+1)*p;  
1330.             F+=p;  
1331.             out++;  
1332.          }  
1333.   
1334.         conclusion.emplace_back(out);  
1335.         i++;  
1336.     }  
1337.     //for_each(conclusion.begin(),conclusion.end(),[](const double&s){cout<<s<<” ”;});  
1338.     //cout<<endl;  
1339.     return c(conclusion);  
1340. }  
1341. //利用间距生成向量  
1342. c seq(const double&from,const double&to,const double&by)  
1343. {  
1344.     if(from>to)  
1345.         throw runtime_error(“from > to\n”);  
1346.     vector<double>out;  
1347.     auto begin=from;  
1348.     while(begin<=to)  
1349.     {  
1350.         out.push_back(begin);  
1351.         begin+=by;  
1352.     }  
1353.     return c(out);  
1354. }  
1355. //生成列表向量元素的所有组合并将结果以矩阵形式输出  
1356. //前面的generate_grid_ 函数为铺垫而设 expand_grid函数接受c为参数给出输出  
1357. matrix generate_grid_0(pair<double,c>p0)  
1358. {  
1359.     matrix out(p0.second.re_dim(),2,0);  
1360.     for(size_t i=0;i<out.re_dim().first;i++)  
1361.     {  
1362.         out[i][0]=p0.first;  
1363.         out[i][1]=p0.second[i];  
1364.     }  
1365.     return out;  
1366. }  
1367. matrix generate_grid_1(const c&c0,const c&c1)  
1368. {  
1369.     vector<double> temp0=c0.const_re_c_val();  
1370.     size_t len=temp0.size();  
1371.     matrix out(len*c1.re_dim(),2,0);  
1372.     for(size_t i=0;i<len;i++)  
1373.     {  
1374.         matrix temp00=generate_grid_0(make_pair(temp0[i],c1));  
1375.        for(size_t j=0;j<c1.re_dim();j++)  
1376.        {  
1377.          out[i*c1.re_dim()+j]=temp00[j];  
1378.        }  
1379.     }  
1380.     return out;  
1381. }  
1382. matrix generate_grid_2(const matrix&m,const c&v)  
1383. {  
1384.     matrix temp=generate_grid_1(m.col(m.re_dim().second-1),v);  
1385.     //temp.out();  
1386.     //system(“pause”);  
1387.     matrix out(temp.re_dim().first,m.re_dim().second+1,0);  
1388.     //out.out();  
1389.     //system(“pause”);  
1390.     matrix add(m.re_dim().first,m.re_dim().second-1,0);  
1391.     for(size_t i=0;i<add.re_dim().second;i++)  
1392.         add.col_change(i,m.col(i));  
1393.     //add.out();  
1394.     //system(“pause”);  
1395.   
1396.     for(size_t i=0;i<m.re_dim().first;i++)  
1397.         for(size_t j=0;j<v.re_dim();j++)  
1398.         {  
1399.             for(size_t k=0;k<add.re_dim().second;k++)  
1400.              out[i*v.re_dim()+j][k]=add[i][k];  
1401.         }  
1402.     out.col_change(out.re_dim().second-2,temp.col(0));  
1403.     out.col_change(out.re_dim().second-1,temp.col(1));  
1404.   
1405.     return out;  
1406. }  
1407. //若干向量其元素的所有组合  
1408. matrix expand_grid(initializer_list<c>i0)  
1409. {  
1410.     vector<c>temp(i0);  
1411.     if(temp.size()==1)  
1412.         return matrix(temp[0]);  
1413.     matrix beg0=generate_grid_1(temp[0],temp[1]);  
1414.     if(temp.size()==2)  
1415.         return beg0;  
1416.     for(size_t i=2;i<temp.size();i++)  
1417.     {  
1418.         //beg0.out();  
1419.         beg0=generate_grid_2(beg0,temp[i]);  
1420.         //cout<<”Done!\n”;  
1421.         //beg0.out();  
1422.         //cout<<”Done!\n”;  
1423.     }  
1424.   
1425.     return beg0;  
1426. }  
1427. //得到向量v的所有n个元素的组合 输出以矩阵排列 每列为一个组合  
1428. //调用了求得下标组合的 require_dim()  
1429. //当组合项数太多时主张利用转置按行查看  
1430. matrix combn(const c&v,size_t n)  
1431. {  
1432.    auto a=require_dim(v.re_dim(),n);  
1433.    auto b=vector<set<size_t>>(a.size());  
1434.    auto a0=a.begin();  
1435.    auto b0=b.begin();  
1436.    while(a0!=a.end())  
1437.    {  
1438.        for_each(a0->begin(),a0->end(),  
1439.                 [&](const pair<size_t,size_t>&p)  
1440.                 {  
1441.                     if(p.second==1)  
1442.                     b0->insert(p.first);  
1443.                 });  
1444.                 a0++;  
1445.                 b0++;  
1446.    }  
1447.     auto cc=vector<vector<double>>(b.size());  
1448.     for(size_t i=0;i<cc.size();i++)  
1449.     {  
1450.         cc[i]=vector<double>(b[i].begin(),b[i].end());  
1451.     }  
1452.       matrix out(cc[0].size(),cc.size(),0);  
1453.       for(size_t i=0;i<out.re_dim().second;i++)  
1454.       {  
1455.          out.col_change(i,cc[i]);  
1456.       }  
1457.       matrix oout(out.re_dim().first,out.re_dim().second,0);  
1458.       for(size_t i=0;i<oout.re_dim().first;i++)  
1459.         for(size_t j=0;j<oout.re_dim().second;j++)  
1460.         oout[i][j]=v[out[i][j]];  
1461.       return oout;  
1462. }  
1463. //上述函数的重载版本 返回0:m的相应n组合  
1464. matrix combn(size_t m,size_t n)  
1465. {  
1466.     vector<double>temp;  
1467.     for(size_t i=0;i<m+1;i++)  
1468.         temp.push_back(i);  
1469.     return combn(temp,n);  
1470. }  
1471. //向量相应元素在前面是否出现  
1472. //出现显示1  
1473. c duplicated(const c&v)  
1474. {  
1475.     //函数a返回向量cc的前n个元素构成的向量  
1476.     auto a=[](const c&cc,size_t n)->c  
1477.     {  
1478.         c out(rep(0,n));  
1479.         for(size_t i=0;i<n;i++)  
1480.             out[i]=cc[i];  
1481.         return out;  
1482.     };  
1483.     vector<c>vc0(v.re_dim());  
1484.     for(size_t i=0;i<vc0.size();i++)  
1485.         vc0[i]=a(v,i+1);  
1486.     vector<vector<double>>vs0(vc0.size());  
1487.     for(size_t i=0;i<vs0.size();i++)  
1488.         vs0[i]=vector<double>(vc0[i].const_re_c_val());  
1489.     c out(rep(0,v.re_dim()));  
1490.     for(size_t i=1;i<v.re_dim();i++)  
1491.     {  
1492.          if(find(vs0[i-1].begin(),vs0[i-1].end(),v[i])!=vs0[i-1].end())  
1493.          out[i]=1;  
1494.     }  
1495.     return out;  
1496. }  
1497. //实现unique的操作看似与factor的levels相同  
1498. //实际在序上是不同的  
1499. //levels保持double的序（由set有序决定）  
1500. //unique应保持原向量的序  
1501. //unique由duplicated定义  
1502. c unique(const c&v)  
1503. {  
1504.     c d0=duplicated(v);  
1505.     vector<pair<double,double>>vm0;  
1506.     for(size_t i=0;i<v.re_dim();i++)  
1507.         vm0.push_back(make_pair(v[i],d0[i]));  
1508.     vector<double>out;  
1509.     for(size_t i=0;i<vm0.size();i++)  
1510.     {  
1511.         if(vm0[i].second==0)  
1512.             out.push_back(vm0[i].first);  
1513.     }  
1514.     return c(out);  
1515. }  
1516. //矩阵叉乘 T(a)%b (向量内积的矩阵版本)  
1517. matrix crossprod(const matrix&a,const matrix&b)  
1518. {  
1519.     return T(a)%b;  
1520. }  
1521. //outer 生成外积 相应的可将外机作用于函数  
1522. matrix outer(const c&v1,const c&v2,function<double(const double&,const double&)>p=function<double(const double&,const double&)>())  
1523. {  
1524.     if(!p)  
1525.     {  
1526.         matrix out(v1.re_dim(),v2.re_dim(),0);  
1527.         for(size_t i=0;i<out.re_dim().first;i++)  
1528.             for(size_t j=0;j<out.re_dim().second;j++)  
1529.             out[i][j]=v1[i]*v2[j];  
1530.         return out;  
1531.     }  
1532.     else  
1533.     {  
1534.         matrix out(v1.re_dim(),v2.re_dim(),0);  
1535.         for(size_t i=0;i<out.re_dim().first;i++)  
1536.             for(size_t j=0;j<out.re_dim().second;j++)  
1537.             out[i][j]=p(v1[i],v2[j]);  
1538.         return out;  
1539.     }  
1540. }  
1541. //矩阵Kronecker积  
1542. matrix kronecker(const matrix&a,const matrix&b)  
1543. {  
1544.     matrix out(a.re_dim().first*b.re_dim().first,a.re_dim().second*b.re_dim().second,0);  
1545.     vector<vector<matrix>>vvm0(a.re_dim().first,vector<matrix>(a.re_dim().second));  
1546.     for(size_t i=0;i<a.re_dim().first;i++)  
1547.         for(size_t j=0;j<a.re_dim().second;j++)  
1548.     {  
1549.        vvm0[i][j]=a[i][j]*b;  
1550.     }  
1551.     for(size_t i=0;i<vvm0.size();i++)  
1552.         for(size_t j=0;j<vvm0[0].size();j++)  
1553.     {  
1554.         matrix temp=vvm0[i][j];  
1555.         for(size_t k=0;k<temp.re_dim().first;k++)  
1556.             for(size_t l=0;l<temp.re_dim().second;l++)  
1557.             out[i*b.re_dim().first+k][j*b.re_dim().second+l]=temp[k][l];  
1558.     }  
1559.     return out;  
1560. }  
1561. //由方阵上下三角部分组成的0-1阵 不包含对角线元素  
1562. matrix lower_tri(const matrix&m)  
1563. {  
1564.     if(m.re_dim().first!=m.re_dim().second)  
1565.         throw runtime_error(“dim error!\n”);  
1566.     matrix out(m.re_dim().first,m.re_dim().second,0);  
1567.     for(size_t i=1;i<out.re_dim().first;i++)  
1568.         for(size_t j=0;j<i;j++)  
1569.         out[i][j]=1;  
1570.     return out;  
1571. }  
1572. matrix upper_tri(const matrix&m)  
1573. {  
1574.     return T(lower_tri(m));  
1575. }  
1576. //矩阵拉直(按列)  
1577. c vec(const matrix&m)  
1578. {  
1579.     c out=m.col(0);  
1580.     for(size_t i=1;i<m.re_dim().second;i++)  
1581.         out=connect({out,m.col(i)});  
1582.     return out;  
1583. }  
1584. //矩阵半拉直 拉直运算忽略主对角线以上部分  
1585. c vech(const matrix&m)  
1586. {  
1587.     c out=m.col(0);  
1588.     for(size_t i=1;i<m.re_dim().second;i++)  
1589.     {  
1590.         vector<double> temp=m.col(i).const_re_c_val();  
1591.         auto a=temp.begin();  
1592.         for(size_t j=0;j<i;j++)  
1593.             a++;  
1594.         vector<double>temp0(a,temp.end());  
1595.         out=connect({out,c(temp0)});  
1596.     }  
1597.     return out;  
1598. }  
1599. //在方阵满秩情况下 给出方阵A的QR分解  
1600. pair<matrix,matrix> qr_of_full_rank(const matrix&A)  
1601. {  
1602.     if(A.re_dim().first!=A.re_dim().second)  
1603.         throw runtime_error(“dim error!\n”);  
1604.     // /*  
1605.     size_t rk=rank(A);  
1606.     cout<<”The rank: ”<<rk<<endl;  
1607.   
1608.     if(rk!=A.re_dim().first)  
1609.         throw runtime_error(“singular!\n”);  
1610.     //*/  
1611.     cout<<”A=QR”<<endl;  
1612.     cout<<endl;  
1613.     cout<<”Q of A:\n”;  
1614.     matrix Q=orth_normal(A);  
1615.     Q.out();  
1616.     cout<<”R of m0:\n”;  
1617.     matrix R=T(Q)%A;  
1618.     R.val_to_double(lower_tri(R),0);  
1619.     R.out();  
1620.     // /*  
1621.     cout<<”The F_norm gap between A and QR:\n”;  
1622.     cout<<F_norm(A-Q%R)<<endl;  
1623.     //*/  
1624.     return make_pair(Q,R);  
1625. }  
1626. //要求自变量矩阵为可逆方阵 利用qr分解 求解线性方程组  
1627. //这只是一个形式上的实现 因为求解R逆的方法是用伴随矩阵法  
1628. //如能改为初等行变换法应会提高效率  
1629. c qr_solve(const matrix&m,const c&v)  
1630. {  
1631.     auto qr=qr_of_full_rank(m);  
1632.     return solve(qr.second)%T(qr.first)%v;  
1633. }  
1634.   
1635. int main()  
1636. {  
1637. //测试程序：  
1638.   
1639.    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);  
1640.    m2.out();  
1641.    c v1({1,1,1,1,2});  
1642.    c v2=solve(m2,v1);  
1643.    c v3=qr_solve(m2,v1);  
1644.    v2.out();  
1645.    v3.out();  
1646.    cout<<”The gap of F_norm:\n”;  
1647.    cout<<F_norm(v2-v3)<<endl;  
1648.    return 0;  
1649. }  
1650.   
1651.    

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;this->operator=this->operator;this->operator=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;
        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法则解线性方程组：要求自变量矩阵为方阵且相应行列式不为0
c 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);
}
//向量相应元素在前面是否出现
//出现显示1
c 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;
}

 

 

转载自斯温的博客 http://blog.csdn.net/sinat_30665603