Advertisement
Draxion

Matrix Operations

Jan 11th, 2017
181
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 35.37 KB | None | 0 0
  1. #include <cstdlib>
  2. #include <iostream>
  3. #include <stdio.h>
  4. #include <ctype.h>
  5. #include <array>
  6. #include <string>
  7. #include <fstream>
  8. #include <sstream>
  9. #include <cmath>
  10. #include "apmatrix.h"
  11.  
  12. void printUM(apmatrix<double> unmod);
  13. void printM(apmatrix<double> matrix);
  14.  
  15. apmatrix<double>add(apmatrix<double>mtx, apmatrix<double>matrix);
  16. apmatrix<double>subtract(apmatrix<double>mtx, apmatrix<double> matrix);
  17. apmatrix<double>mult(apmatrix<double>mtx, apmatrix<double> matrix);
  18. apmatrix<double>rE(apmatrix<double> matrix, apmatrix<double>unmod);
  19. apmatrix<double>rrE(apmatrix<double> matrix, apmatrix<double>unmod);
  20. apmatrix<double>transpose(apmatrix<double>matrix, apmatrix<double>unmod);
  21. apmatrix<double>inverse(apmatrix<double> mtx, apmatrix<double>unmod);
  22. apmatrix<double>linReg(apmatrix<double>x,apmatrix<double>y);
  23. apmatrix<double>regrssn(apmatrix<double>x,apmatrix<double>y, apmatrix<string> f, apmatrix<double> v);
  24. apmatrix<double>augment(apmatrix<double>A,apmatrix<double>B);
  25. apmatrix<double>redaug(apmatrix<double>aug, apmatrix<double> A, apmatrix<double>B);
  26. apmatrix<double>proj(apmatrix<double>v, apmatrix<double>u);
  27. apmatrix<double>orthogonal(apmatrix<double>A);
  28. apmatrix<double>orthonormal(apmatrix<double>V);
  29. void enter();
  30. void setMtx(int indx, apmatrix<double>mtx);
  31. void store(apmatrix<double>matrix,apmatrix<double>orig);
  32. void rVals();
  33. apmatrix<double> getMtx(int val);
  34. apmatrix<double> getUnM(int val);
  35. apmatrix<double> matrix,m1,m2,m3,m4,m5,m6,unmod,unmod1, unmod2,unmod3,unmod4,unmod5,unmod6, dataX, dataY, dataV;
  36. apmatrix<string> dataF;
  37. int row, col;
  38. using namespace std;
  39.  
  40. int main(int argc, char *argv[])
  41. {
  42.     bool menu=true, loop = false, tp = false, red = false;
  43.     while(menu) //Menu
  44.     {
  45.         cout<<"\nMenu\n"<<"----\n";
  46.         cout<<"1) Enter Matrix\n"<<"2) Show Matrix\n"<<"3) Transpose Matrix\n"<<"4) Add Matrices\n"<<"5) Subtract Matrices\n"<<"6) Multiply Matrices\n"<<"7) Row Echelon\n"<<"8) Reduced Row Echelon\n"<<"9) Inverse Matrix\n"<<"10) Linear Regression\n"<<"11) Matrix Regression\n"<<"12) QR Factorization\n"<<"0) Quit\n"<<"----\n>";
  47.         int val;
  48.         cin>>val;
  49.         if(val == 1)
  50.         {
  51.             enter();
  52.             system("PAUSE");
  53.         }
  54.         else if(val==2)   //Show Matrix
  55.         {
  56.             bool l;
  57.             while(l)
  58.             {
  59.                 cout<<"Which matrix do you want to see (1-6)? "<<"\n>";
  60.                 int show;
  61.                 cin>>show;
  62.                 if(show>=1 && show<=6)
  63.                 {
  64.                     printM(getMtx(show));
  65.                     l=false;
  66.                 }
  67.                 else
  68.                 {
  69.                     cout<<"Please enter a positive integer from 1-6.\n";
  70.                 }
  71.             }
  72.             cout<<"\n";
  73.             system("PAUSE");
  74.         }
  75.         else if(val == 3) //Transpose Matrix
  76.         {
  77.             int num;
  78.             cout<<"Which matrix do you want to transpose?"<<"\n>";
  79.             cin>>num;
  80.             //Changes original matrix to matrix stored in num value.
  81.             setMtx(num,getMtx(num));
  82.             matrix = transpose(getMtx(num), getUnM(num));
  83.             printUM(getUnM(num));
  84.             cout<<"Transposed Matrix:"<<"\n";
  85.             printM(matrix);
  86.             store(matrix, getUnM(num));
  87.             system("PAUSE");
  88.         }
  89.         else if(val == 4) //Add Matrices
  90.         {
  91.             int mtx1, mtx2;
  92.             apmatrix<double>sum;
  93.             cout<<"Enter first matrix to add: ";
  94.             cin>>mtx1;
  95.             cout<<"Enter second matrix to add: ";
  96.             cin>>mtx2;
  97.             sum = add(getMtx(mtx1), getMtx(mtx2));
  98.             if(sum[1][1] == 1)
  99.             {
  100.                 cout<<"Matrix Size Mismatch. <"<<getMtx(mtx1).numrows()<<","<<getMtx(mtx1).numcols()<<"> + <"<<getMtx(mtx2).numrows()<<","<<getMtx(mtx2).numcols()<<"> is not possible.\n";
  101.                 system("PAUSE");
  102.             }
  103.             else
  104.             {
  105.                 printM(getMtx(mtx1));
  106.                 //cout<<"\n";
  107.                 for(int i=0; i<getMtx(mtx1).numcols()/2; i++)
  108.                     cout<<" ";
  109.                 cout<<"+\n";
  110.                 printM(getMtx(mtx2));
  111.                 cout<<"\n-";
  112.                 for(int i=0; i<sum.numcols(); i++)
  113.                     cout<<"-";
  114.                 cout<<"\n";
  115.                 printM(sum);
  116.                 store(sum, sum);
  117.             }
  118.             system("PAUSE");
  119.         }
  120.         else if(val == 5) //Subtract Matrices
  121.         {
  122.             int mtx1, mtx2;
  123.             apmatrix<double>diff;
  124.             cout<<"Enter first matrix to subtract: ";
  125.             cin>>mtx1;
  126.             cout<<"Enter second matrix to subtract: ";
  127.             cin>>mtx2;
  128.             diff = subtract(getMtx(mtx1), getMtx(mtx2));
  129.             if(diff[1][1] == 1)
  130.             {
  131.                 cout<<"Matrix Size Mismatch. <"<<getMtx(mtx1).numrows()<<","<<getMtx(mtx1).numcols()<<"> - <"<<getMtx(mtx2).numrows()<<","<<getMtx(mtx2).numcols()<<"> is not possible.\n";
  132.                 system("PAUSE");
  133.             }
  134.             else
  135.             {
  136.                 printM(getMtx(mtx1));
  137.                 //cout<<"\n";
  138.                 for(int i=0; i<getMtx(mtx1).numcols()/2; i++)
  139.                     cout<<" ";
  140.                 cout<<"-\n";
  141.                 printM(getMtx(mtx2));
  142.                 cout<<"\n-";
  143.                 for(int i=0; i<diff.numcols(); i++)
  144.                     cout<<"-";
  145.                 cout<<"\n";
  146.                 printM(diff);
  147.                 store(diff, diff);
  148.             }
  149.             system("PAUSE");
  150.         }
  151.         else if(val == 6) //Multiply Matrices
  152.         {
  153.             int mtx1, mtx2;
  154.             apmatrix<double>prod;
  155.             cout<<"Enter first matrix to multiply: ";
  156.             cin>>mtx1;
  157.             cout<<"Enter second matrix to multiply: ";
  158.             cin>>mtx2;
  159.             prod = mult(getMtx(mtx1), getMtx(mtx2));
  160.             if(prod[1][1] == 1)
  161.             {
  162.                 cout<<"Matrix Size Mismatch. <"<<getMtx(mtx1).numrows()<<","<<getMtx(mtx1).numcols()<<"> x <"<<getMtx(mtx2).numrows()<<","<<getMtx(mtx2).numcols()<<"> is not possible.\n";
  163.                 system("PAUSE");
  164.             }
  165.             else
  166.             {
  167.                 printM(getMtx(mtx1));
  168.                 //cout<<"\n";
  169.                 for(int i=0; i<getMtx(mtx1).numcols()/2; i++)
  170.                     cout<<" ";
  171.                 cout<<"*\n";
  172.                 printM(getMtx(mtx2));
  173.                 cout<<"\n-";
  174.                 for(int i=0; i<prod.numcols(); i++)
  175.                     cout<<"-";
  176.                 cout<<"\n";
  177.                 printM(prod);
  178.                 store(prod, prod);
  179.             }
  180.             system("PAUSE");
  181.         }
  182.         else if(val == 7) //Row Echelon
  183.         {
  184.             int num;
  185.             cout<<"Which matrix do you want to reduce?\n>";
  186.             cin>>num;
  187.             //Changes original/unmod matrix to matrix stored in num value.
  188.             setMtx(num,getMtx(num));
  189.             printUM(getUnM(num));
  190.             matrix.resize(getMtx(num).numrows(), getMtx(num).numcols()); //resize to fit target matrix
  191.             matrix = rE(getMtx(num), getUnM(num)); //matrix = target matrix
  192.             cout<<"REF Matrix:"<<"\n";
  193.             printM(matrix); //Print new matrix
  194.             store(matrix,getUnM(num)); //Store matrix
  195.             system("PAUSE");
  196.         }
  197.         else if(val == 8) //Reduced Row Echelon
  198.         {
  199.             int num;
  200.             cout<<"Which matrix do you want to reduce?\n>";
  201.             cin>>num;
  202.             //Changes original/unmod matrix to matrix stored in num value.
  203.             setMtx(num,getMtx(num));
  204.             printUM(getUnM(num));
  205.             matrix.resize(getMtx(num).numrows(), getMtx(num).numcols()); //resize to fit target matrix
  206.             matrix = rrE(rE(getMtx(num), getUnM(num)),rE(getMtx(num), getUnM(num))); //matrix = target matrix
  207.             cout<<"RREF Matrix:"<<"\n";
  208.             printM(matrix); //Print new matrix
  209.             store(matrix,getUnM(num)); //Store matrix
  210.             system("PAUSE");
  211.         }
  212.         else if(val == 9) //Inverse Matrix
  213.         {
  214.             int num;
  215.             cout<<"Which matrix do you want to invert?\n>";
  216.             cin>>num;
  217.             //Changes original/unmod matrix to matrix stored in num value.
  218.             setMtx(num,getMtx(num));
  219.             printUM(getUnM(num));
  220.             matrix.resize(getMtx(num).numrows(), getMtx(num).numcols()); //resize to fit target matrix
  221.             matrix = inverse(getMtx(num),getUnM(num)); //matrix = target matrix
  222.             cout<<"Inverse Matrix:"<<"\n";
  223.             printM(matrix); //Print new matrix
  224.             store(matrix,getUnM(num)); //Store matrix
  225.             system("PAUSE");
  226.         }
  227.         else if(val == 10) //Linear Regression
  228.         {
  229.             rVals();
  230.             //printM(dataY);
  231.             apmatrix<double>reg;
  232.             reg.resize(1,2);
  233.             reg = linReg(dataX, dataY);
  234.             cout<<"y=ax+b => y="<<reg[0][1]<<"x+"<<reg[0][0];
  235.             system("PAUSE");
  236.         }
  237.         else if(val == 11) //Matrix Regression
  238.         {
  239.             rVals();
  240.             cout<<"Enter operations. Type 'm' for scalar multiplication, 'p' for power, 'sqrt' for square root, 'l' for log, 'sin',  'cos', 'tan', or 'd' when finished.\n";
  241.             bool eF = true;
  242.             int row=0;
  243.             dataF.resize(row+1, 1);
  244.             dataV.resize(row+1, 1);
  245.             while(eF)
  246.             {
  247.                 cout<<">";
  248.                 string n;
  249.                 cin>>n;
  250.                 if(n=="d")
  251.                 {
  252.                     dataF.resize(dataV.numrows()-1,1);
  253.                     dataV.resize(dataF.numrows(),1);
  254.                     eF = false;
  255.                 }
  256.  
  257.                 else
  258.                 {
  259.                     if(n=="m")
  260.                     {
  261.                         cout<<"Enter scalar value: ";
  262.                         cin>>dataV[row][0];
  263.                     }
  264.                     else if(n == "p")
  265.                     {
  266.                         cout<<"Enter exponent: ";
  267.                         cin>>dataV[row][0];
  268.                     }
  269.                     else if(n == "l")
  270.                     {
  271.                         cout<<"Enter log base value: ";
  272.                         cin>>dataV[row][0];
  273.                     }
  274.                     else
  275.                     {
  276.                         dataV[row][0] = 0;
  277.                     }
  278.                     double num;
  279.                     // dataF[row][0]=strtod(n.c_str(),NULL);
  280.                     dataF[row][0] = n;
  281.                     row++;
  282.                     dataF.resize(row+1,1);
  283.                     dataV.resize(row+1,1);
  284.                 }
  285.             }
  286.             apmatrix<double> view = regrssn(dataX, dataY, dataF, dataV);
  287.             /*for(int i = 0; i<view.numrows(); i++) {
  288.                 for(int j = 0; j<view.numcols(); j++) {
  289.                     cout<<view[i][j];
  290.                 }
  291.                 cout<<endl;
  292.             }*/
  293.             cout<<"y="<<view[0][0]<<"+"<<view[1][0]<<"x+";
  294.             for(int i = 2; i<view.numrows(); i++)
  295.             {
  296.                 for(int j = 0; j<view.numcols(); j++)
  297.                 {
  298.                     if(dataF[i-2][0]==("m") || dataF[i-2][0]==("M"))   //Scalar multiplication
  299.                     {
  300.                         cout<<view[i][0]<<"*x";
  301.                     }
  302.                     else if(dataF[i-2][0]==("p") || dataF[i-2][0]==("P"))   //Power
  303.                     {
  304.                         cout<<view[i][0]<<"*";
  305.                         cout<<"x^"<<dataV[i-2][0];
  306.                     }
  307.                     else if(dataF[i-2][0]==("sqrt"))   //Square Root
  308.                     {
  309.                         cout<<view[i][0]<<"*sqrt(x)";
  310.                     }
  311.                     else if (dataF[i-2][0]==("l") || dataF[i-2][0]==("L"))   //Log Base
  312.                     {
  313.                         cout<<view[i][0]<<"*log("<<dataV[i-2][0]<<",x)";
  314.                     }
  315.                     else if (dataF[i-2][0]==("sin"))   //Sin
  316.                     {
  317.                         cout<<view[i][0]<<"*sin(x)";
  318.                     }
  319.                     else if (dataF[i-2][0]==("cos"))   //Cos
  320.                     {
  321.                         cout<<view[i][0]<<"*cos(x)";
  322.                     }
  323.                     else if (dataF[i-2][0]==("tan"))   //Tan
  324.                     {
  325.                         cout<<view[i][0]<<"*tan(x)";
  326.                     }
  327.                 }
  328.                 if(i+2<view.numrows()&&view[i+1][0]>=0)
  329.                 {
  330.                     cout<<"+";
  331.                 }
  332.             }
  333.         }
  334.         else if(val == 12) //QR Factorization
  335.         {
  336.             int num;
  337.             cout<<"Which matrix do you want to factorize?\n>";
  338.             cin>>num;
  339.             //Changes original/unmod matrix to matrix stored in num value.
  340.             setMtx(num,getMtx(num));
  341.             printUM(getUnM(num));
  342.             apmatrix<double>Q = orthonormal(orthogonal(getUnM(num)));
  343.             apmatrix<double>R = mult(transpose(Q,Q),getUnM(num));
  344.             cout<<"Q:\n";
  345.             for(int i = 0; i<Q.numrows(); i++)
  346.             {
  347.                 for(int j=0; j<Q.numcols(); j++)
  348.                 {
  349.                     cout<<Q[i][j]<<" ";
  350.                 }
  351.                 cout<<endl;
  352.             }
  353.             cout<<"\nR:\n";
  354.             for(int i = 0; i<R.numrows(); i++)
  355.             {
  356.                 for(int j=0; j<R.numcols(); j++)
  357.                 {
  358.                     cout<<R[i][j]<<" ";
  359.                 }
  360.                 cout<<endl;
  361.             }
  362.             system("PAUSE");
  363.         }
  364.         else if(val == 0) //Quit
  365.         {
  366.             // onClose();
  367.             return EXIT_SUCCESS;
  368.         }
  369.         else
  370.         {
  371.             cout<<"Please enter a 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, or 0.\n";
  372.         }
  373.     }//I
  374.     system("PAUSE");
  375.     return EXIT_SUCCESS;
  376. }
  377. void enter()
  378. {
  379.     //Declare Matrix Size
  380.     apmatrix<double>mtx;
  381.     apmatrix<double>orig;
  382.     cout<<"Enter the 'm' size of the matrix: "<<"\t";
  383.     cin>>row;
  384.     cout<<"Enter the 'n' size of the matrix: "<<"\t";
  385.     cin>>col;
  386.     //Enter Values for Matrix
  387.     mtx.resize(row, col);
  388.     orig.resize(row,col);
  389.     for(int f=0; f<row; f++) //A
  390.     {
  391.         for(int j=0; j<col; j++) //1
  392.         {
  393.             double coord;
  394.             cout<<"Enter value for ("<<(f+1)<<","<<(j+1)<<"): ";
  395.             cin>>coord;
  396.             mtx[f][j]=coord;
  397.             orig[f][j]=coord;
  398.         } //1
  399.     } //A
  400.     //Overwrite or Store
  401.     bool st;
  402.     while(st)
  403.     {
  404.         cout<<"Which matrix would you like to store this in (1-6)? Enter '0' to go to menu.\n>";
  405.         int indx;
  406.         cin>>indx;
  407.         if(indx>=1 && indx<=6)
  408.         {
  409.             setMtx(indx,mtx);
  410.             cout<<"\nThis matrix is labeled: "<<indx<<"\n";
  411.  
  412.             st = false;
  413.         }
  414.         else if(indx==0)
  415.             st=false;
  416.         else
  417.         {
  418.             cout<<"\nPlease enter a whole number value from 1-6\n";
  419.         }
  420.     }
  421. }
  422. void setMtx(int indx, apmatrix<double>mtx)
  423. {
  424.     if(indx == 1)
  425.     {
  426.         m1.resize(row,col);
  427.         m1 = mtx;
  428.         unmod1.resize(row,col);
  429.         unmod1 = mtx;
  430.     }
  431.     else if(indx == 2)
  432.     {
  433.         m2.resize(row,col);
  434.         m2=mtx;
  435.         unmod2.resize(row,col);
  436.         unmod2 = mtx;
  437.     }
  438.     else if(indx == 3)
  439.     {
  440.         m3.resize(row,col);
  441.         m3 = mtx;
  442.         unmod3.resize(row,col);
  443.         unmod3 = mtx;
  444.     }
  445.     else if(indx == 4)
  446.     {
  447.         m4.resize(row,col);
  448.         m4 = mtx;
  449.         unmod4.resize(row,col);
  450.         unmod4 = mtx;
  451.     }
  452.     else if(indx == 5)
  453.     {
  454.         m5.resize(row,col);
  455.         m5 = mtx;
  456.         unmod5.resize(row,col);
  457.         unmod5 = mtx;
  458.     }
  459.     else if (indx== 6)
  460.     {
  461.         m6.resize(row,col);
  462.         m6 = mtx;
  463.         unmod6.resize(row,col);
  464.         unmod6 = mtx;
  465.     }
  466.     else
  467.     {
  468.         cout<<"Something is not right.";
  469.         system("EXIT_FAILURE");
  470.     }
  471. }
  472. apmatrix<double> getMtx(int val)
  473. {
  474.     if (val ==1)
  475.     {
  476.         return m1;
  477.     }
  478.     else if(val ==2)
  479.     {
  480.         return m2;
  481.     }
  482.     else if (val ==3)
  483.     {
  484.         return m3;
  485.     }
  486.     else if(val==4)
  487.     {
  488.         return m4;
  489.     }
  490.     else if(val==5)
  491.     {
  492.         return m5;
  493.     }
  494.     else if(val=6)
  495.     {
  496.         return m6;
  497.     }
  498. }
  499. apmatrix<double> getUnM(int val)
  500. {
  501.     if (val ==1)
  502.     {
  503.         return unmod1;
  504.     }
  505.     else if(val ==2)
  506.     {
  507.         return unmod2;
  508.     }
  509.     else if (val ==3)
  510.     {
  511.         return unmod3;
  512.     }
  513.     else if(val==4)
  514.     {
  515.         return unmod4;
  516.     }
  517.     else if(val==5)
  518.     {
  519.         return unmod5;
  520.     }
  521.     else if(val=6)
  522.     {
  523.         return unmod6;
  524.     }
  525. }
  526. void store(apmatrix<double>matrix, apmatrix<double>orig)
  527. {
  528.     bool h =true;
  529.     while(h)
  530.     {
  531.         cout<<"Which matrix would you like to store this in?\nEnter 'c' for current, '0' to go to menu, or enter 1-6.\n>";
  532.         char stor;
  533.         cin>>stor;
  534.         if(stor == 'c' || stor=='C')   //Stores matrix in state that it's already in.
  535.         {
  536.             bool test = true;
  537.             for(int i=1; i<=6; i++) //Iterates through each matrix
  538.             {
  539.                 if(getMtx(i).numrows() == orig.numrows() && getMtx(i).numcols() == orig.numcols())   //Checks that dimensions are equal
  540.                 {
  541.                     for(int k = 0; k<orig.numrows(); k++)  //Unable to test condition: if(getMatrix(i) == orig)
  542.                     {
  543.                         for(int j = 0; j<orig.numcols(); j++)
  544.                         {
  545.                             if(getMtx(i)[k][j]!=orig[k][j]) //Individually checks each index until one does not match
  546.                                 test=false;
  547.                         }
  548.                     }
  549.                 }
  550.                 if(test)   //If all indexes are same, overwrites matrix and exits loop
  551.                 {
  552.                     setMtx(i,matrix);
  553.                     h=false;
  554.                 }
  555.             }
  556.             if(!test)
  557.             {
  558.                 cout<<"[WARNING] This matrix does not exist in any existing slot. Please choose another option.";
  559.             }
  560.         }
  561.         else if(stor>=49 && stor<=54)     //Int value of the char '1' - '6' (Look at ASCII table)
  562.         {
  563.             setMtx(stor-49,matrix);
  564.             h=false;
  565.         }
  566.         else if(stor=='0') //quit
  567.         {
  568.             h=false;
  569.         }
  570.         else
  571.         {
  572.             cout<<"Please enter 'c' for current, '0' to go to menu or enter 1-6\n";
  573.         }
  574.     }
  575. }
  576.  
  577. apmatrix<double> transpose(apmatrix<double> matrix, apmatrix<double>orig)    //orig = unmod matrix
  578. {
  579.     apmatrix<double>trans;
  580.     trans.resize(matrix.numcols(), matrix.numrows());
  581.     // matrix.resize(orig.numrows(),orig.numcols()); //If not a square, flips m with n
  582.     for(int i = 0; i<matrix.numrows(); i++)
  583.     {
  584.         for(int j=0; j<matrix.numcols(); j++)
  585.         {
  586.  
  587.             trans[j][i] = matrix[i][j]; //Switching the variables
  588.         }
  589.     }
  590.     return trans;
  591. }
  592. apmatrix<double> add(apmatrix<double>mtx, apmatrix<double> matrix)
  593. {
  594.     apmatrix<double> ph;
  595.     if(mtx.numrows() == matrix.numrows() && mtx.numcols() == matrix.numcols())
  596.     {
  597.         for(int i = 0; i<matrix.numrows(); i++)
  598.         {
  599.             for(int j = 0; j<matrix.numcols(); j++)
  600.             {
  601.                 ph.resize(matrix.numrows(), matrix.numcols());
  602.                 ph[i][j] = mtx[i][j] + matrix[i][j];
  603.             }
  604.         }
  605.         return ph;
  606.  
  607.     }
  608.     else
  609.     {
  610.         apmatrix<double>m;
  611.         // m.resize(1,1);
  612.         m[0][0]=1.0;
  613.         return m;
  614.     }
  615.     // return matrix;
  616. }
  617. apmatrix<double> subtract(apmatrix<double>mtx, apmatrix<double> matrix)
  618. {
  619.     apmatrix<double> ph;
  620.     if(mtx.numrows() == matrix.numrows() && mtx.numcols() == matrix.numcols())
  621.     {
  622.         for(int i = 0; i<matrix.numrows(); i++)
  623.         {
  624.             for(int j = 0; j<matrix.numcols(); j++)
  625.             {
  626.                 ph.resize(matrix.numrows(), matrix.numcols());
  627.                 ph[i][j] = mtx[i][j] - matrix[i][j];
  628.             }
  629.         }
  630.         return ph;
  631.  
  632.     }
  633.     else
  634.     {
  635.         apmatrix<double>m;
  636.         // m.resize(1,1);
  637.         m[0][0]=1.0;
  638.         return m;
  639.     }
  640.     // return matrix;
  641. }
  642. apmatrix<double> mult(apmatrix<double>mtx, apmatrix<double> matrix)
  643. {
  644.     apmatrix<double>dh;
  645.     if(mtx.numcols() == matrix.numrows())   //Checks if col of mtx1 = row of mtx2
  646.     {
  647.         dh.resize(mtx.numrows(), matrix.numcols()); //Sets size (mtx1.col, mtx2.row);
  648.         // cout<<"rows:rows: "<<dh.numrows()<<" cols: "<<dh.numcols()<<"\n";
  649.         for(int i = 0; i<dh.numrows(); i++)
  650.         {
  651.             for(int j = 0; j<dh.numcols(); j++)   //This loop is for value to be entered
  652.             {
  653.                 //double index=0;
  654.                 for(int k=0; k<mtx.numcols(); k++)
  655.                 {
  656.                     dh[i][j]+=mtx[i][k]*matrix[k][j];
  657.                 }
  658.                 //dh[i][j] = index;
  659.             }
  660.         }
  661.         return dh;
  662.  
  663.     }
  664.     else
  665.     {
  666.         apmatrix<double>m;
  667.         // m.resize(1,1);
  668.         m[0][0]=1.0;
  669.         return m;
  670.     }
  671. }
  672. apmatrix<double> rE(apmatrix<double> matrix, apmatrix<double> unmd)
  673. {
  674. //Row Echelon Form
  675.     /*Variables:*/
  676.     //cout<<"Matrix after import: "; printM(matrix);
  677.     for(int r=0; r<matrix.numrows()-1; r++) //B //Target Row
  678.     {
  679.         for(int i =0; i<matrix.numrows(); i++)   //Checks if all leading values are 0.
  680.         {
  681.             if(matrix[i][0]!=0)
  682.             {
  683.                 i=matrix.numrows()-1;
  684.                 continue;
  685.             }
  686.             else if (i==matrix.numrows()-1 && matrix[i][0]==0)
  687.                 r++;
  688.         }
  689.         if(matrix[r][r]==0)  //1 //Leading 0 Checker
  690.         {
  691.             bool t = true;
  692.             while(t)  //a //Instead of break statement.
  693.             {
  694.                 for(int x=r+1; x<matrix.numrows(); x++) //i //Switches Target row if leading coefficient is 0
  695.                 {
  696.                     if(matrix[x][r]!=0)  //! //Main Diagonal
  697.                     {
  698.                         double temp[matrix.numcols()];//=matrix[r];
  699.                         for(int h=0; h<matrix.numcols(); h++) //z //Flips rows if leading value is 0
  700.                         {
  701.                             temp[h] = matrix[0][h];
  702.                             matrix[r][h]=matrix[x][h];
  703.                             matrix[x][h]=temp[h];
  704.                         }//z
  705.                         t=false;
  706.                     }//!
  707.                 }//i
  708.             }//a
  709.         }//1
  710.         // cout<<"Matrix after 0 checker: "; printM(matrix);
  711.         for(int i=r+1; i<matrix.numrows(); i++) //b
  712.         {
  713.             for(int j=matrix.numcols()-1; j>=0; j--) //i //Simplification
  714.             {
  715.                 if(matrix[i][r]!=0)  //!
  716.                 {
  717.                     //Creates leading 0
  718.                     //   cout<<"Math: "<<matrix[i][j]<<"="<<matrix[r][j]<<"*"<<matrix[i][r]<<"/"<<matrix[r][r]<<"-"<<matrix[i][j];
  719.                     matrix[i][j]=(double)((double)(matrix[r][j]*(double)(matrix[i][r]/matrix[r][r]))-(double)matrix[i][j]); //Rt*(Rm/Rt)-Rm
  720.                     // cout<<"\nr: "<<r<<"\ni: "<<i<<"\nj: "<<j;
  721.                     // cout<<"="<<matrix[i][j]<<endl;
  722.                 }//!
  723.  
  724.             } //i
  725.         } //b
  726.  
  727.     } //B
  728.  
  729.     //Makes leading value 1
  730.     for(int i = 0; i<matrix.numrows(); i++)
  731.     {
  732.         for(int j = matrix.numcols()-1; j>=i; j-- )
  733.         {
  734.  
  735.             if(matrix[i][i]!=0)
  736.             {
  737.                 //cout<<"Leading 1: "<<matrix[i][j]<<"="<<matrix[i][j]<<"/"<<matrix[i][i]<<"=";
  738.                 double temp = (double)(matrix[i][j]/matrix[i][i]);
  739.                 matrix[i][j] = temp;
  740.                 //cout<<matrix[i][j]<<endl;
  741.             }
  742.         }
  743.     }
  744.     return matrix;
  745. }
  746. apmatrix<double> rrE(apmatrix<double> matrix, apmatrix<double> unmod)
  747. {
  748.     for(int r=0; r<matrix.numrows()-1; r++) //Targeted Row
  749.     {
  750.         for(int i=r+1; i<matrix.numrows(); i++) //Alter row
  751.         {
  752.             for(int j=r+1; j<matrix.numcols(); j++) //Col of Alter row
  753.             {
  754.                 //if(unmod[i][r+1]!=0) { //Makes sure leading value in targeted !=0
  755.  
  756.                 matrix[r][j] = (matrix[r][j]-unmod[i][j]*unmod[r][r+1]/unmod[i][r+1]); //B=B-E(Rt/Et)
  757.                 //}
  758.             }
  759.         }
  760.     }
  761.     return matrix;
  762. }
  763. apmatrix<double> inverse(apmatrix<double> mtx, apmatrix<double> unmd)
  764. {
  765.     apmatrix<double>invert;
  766.     invert.resize(mtx.numrows(),mtx.numcols());
  767.     mtx.resize(mtx.numrows(),(2*mtx.numcols())); //Space to add Identity
  768. //Adds Identity to end of matrix
  769.     for(int i = 0; i<mtx.numrows(); i++)
  770.     {
  771.         for(int j = (mtx.numcols()/2); j<mtx.numcols(); j++)
  772.         {
  773.             if(i==j-(mtx.numcols()/2))   //Finding where to put leading 1
  774.             {
  775.                 mtx[i][j]=1;
  776.                 //  invert[i][j-mtx.numcols()/2]=1; //Becomes identity
  777.             }
  778.             else
  779.             {
  780.                 mtx[i][j] = 0; //If not spot for leading 1, set index to 0.
  781.                 // invert[i][j-(mtx.numcols()/2)]=0;
  782.             }
  783.         }
  784.     }
  785. //matrix = mtx;
  786. // mtx=rE(mtx,mtx);
  787.     //cout<<"RE: "; printM(rE(mtx,mtx));
  788.     mtx = rrE(rE(mtx,unmd),rE(mtx,unmd));//Reduces to get inverse on other side
  789.     for(int i = 0; i<mtx.numrows(); i++)
  790.     {
  791.         for(int j = (mtx.numcols()/2); j<mtx.numcols(); j++)
  792.         {
  793.             invert[i][j-(mtx.numcols()/2)]=mtx[i][j]; //Sets invert to inverted matrix
  794.         }
  795.     }
  796.  
  797.     return invert;
  798. }
  799. //Prints matrix before operation
  800. void printUM(apmatrix<double> unmod)
  801. {
  802.     int row = unmod.numrows();
  803.     int col = unmod.numcols();
  804.  
  805.     cout<<"\n"<<"Original Matrix:"<<"\n"<<"[";
  806.     //Formatting
  807.     for(int i=0; i<row; i++) //D
  808.     {
  809.         for(int j=0; j<col; j++) //1
  810.         {
  811.             if(j==0&&i!=0)  //a
  812.             {
  813.                 cout<<" ";
  814.             }//a
  815.             cout<<unmod[i][j];
  816.  
  817.             if(j<col-1) //b
  818.             {
  819.                 cout<<" ";
  820.             }//b
  821.             else if (i<row-1)  //c
  822.             {
  823.                 cout<<"\n";
  824.             }//c
  825.         }//2
  826.     }//D
  827.     cout<<"]";
  828.     cout<<"\n"<<"\n";
  829. }
  830. //Prints matrix after operation
  831. void printM(apmatrix<double> matrix)
  832. {
  833.     //Print out REF Matrix
  834.     int row = matrix.numrows();
  835.     int col = matrix.numcols();
  836.  
  837.     cout<<"[";
  838.     for(int i=0; i<row; i++) //E
  839.     {
  840.         for(int j=0; j<col; j++) //1
  841.         {
  842.             if(j==0&&i!=0)  //a
  843.             {
  844.                 cout<<" ";
  845.             }//a
  846.             cout<<matrix[i][j];
  847.  
  848.             if(j<col-1) //b
  849.             {
  850.                 cout<<" ";
  851.             }//b
  852.             else if (i<row-1)  //c
  853.             {
  854.                 cout<<"\n";
  855.             }//c
  856.         }//1
  857.     }//E
  858.     cout<<"]";
  859.     cout<<"\n";
  860. }
  861. apmatrix<double> linReg(apmatrix<double>x, apmatrix<double>y)   //Linear Regression
  862. {
  863.     double a=0,b=0,sigX=0,sigY=0,sigX2=0, sigXY=0;
  864.     int n = x.numrows(); //Sig = sigma
  865.     for(int i = 0; i<x.numrows(); i++)
  866.     {
  867.         for(int j = 0; j<y.numcols(); j++)   //Only 1 column
  868.         {
  869.             sigX+=x[i][j]; //Sum of all x values
  870.             sigX2+=pow(x[i][j],2); //Sum of all x^2 valued.
  871.             sigY+=y[i][j]; //Sum of all y values
  872.             sigXY+=x[i][j]*y[i][j]; //Sum of all x*y values
  873.         }
  874.     }
  875.     b=(double)(sigY*sigX2-sigX*sigXY)/((double)n*sigX2-pow(sigX,2)); //formula for y-intercept
  876.     a=(double)((double)n*sigXY-sigX*sigY)/((double)n*sigX2-pow(sigX,2)); //formula for slope
  877.     apmatrix<double>test;
  878.     test.resize(1,2); //Just for logistic purposes
  879.     test[0][0]=a;
  880.     test[0][1]=b;
  881.     return test; //Inefficient
  882. }
  883. apmatrix<double> regrssn(apmatrix<double> x, apmatrix<double> y, apmatrix<string> f/*function*/, apmatrix<double> v /*degree of f*/)   //Z = [(MtM)^-1] * Mt * Y
  884. {
  885.     apmatrix<double>m; //Gets rid of overwrite bug.
  886.     //x.resize(x.numrows(), x.numcols()+1);
  887.     m.resize(x.numrows(), f.numrows()+2);
  888.  
  889.     for(int i = 0; i<x.numrows(); i++)   //Adds row of 1s to x matrix.
  890.     {
  891.         m[i][0] = 1;
  892.         m[i][1] = x[i][0];
  893.         for(int j = 2; j<m.numcols(); j++)   //Placeholder for uninitialized columns
  894.         {
  895.             m[i][j]=0;
  896.         }
  897.     }
  898.     for(int i = 0; i<f.numrows(); i++)   //Number of functions to add to M matrix
  899.     {
  900.         //   x.resize(x.numrows(), x.numcols()+1); //Adds another column to x functions (M)
  901.         for(int j = 0; j<x.numrows(); j++)
  902.         {
  903.             if(f[i][0]==("m") || f[i][0]==("M"))   //Scalar multiplication
  904.             {
  905.                 m[j][i+2] =  v[i][0]*x[j][0];
  906.             }
  907.             else if(f[i][0]==("p") || f[i][0]==("P"))   //Power
  908.             {
  909.                 m[j][i+2] = pow(x[j][0],v[i][0]);
  910.             }
  911.             else if(f[i][0]==("sqrt"))   //Square Root
  912.             {
  913.                 m[j][i+2] = sqrt(x[j][0]);
  914.             }
  915.             else if (f[i][0]==("l") || f[i][0]==("L"))   //Log Base
  916.             {
  917.                 m[j][i+2] = log(x[j][0])/log(v[i][0]);
  918.             }
  919.             else if (f[i][0]==("sin"))   //Sin
  920.             {
  921.                 m[j][i+2] = sin(x[j][0]);
  922.             }
  923.             else if (f[i][0]==("cos"))   //Cos
  924.             {
  925.                 m[j][i+2] = cos(x[j][0]);
  926.             }
  927.             else if (f[i][0]==("tan"))   //Tan
  928.             {
  929.                 m[j][i+2] = tan(x[j][0]);
  930.             }
  931.             else
  932.             {
  933.                 continue;
  934.             }
  935.         }
  936.     }
  937.     apmatrix<double> M = m;
  938.     apmatrix<double> Mt = transpose(m,m);
  939.     apmatrix<double> MtM = mult(Mt, M);
  940.     y.resize(y.numrows()-1,y.numcols()); //Extra row of a 0 for some reason.
  941.     apmatrix<double>MtY = mult(Mt,y);
  942.     //apmatrix<double> MtMI = inverse(MTM,MTM);
  943.  
  944.     // cout<<"M: "<<M.numrows()<<"x"<<M.numcols()<<endl;
  945.     // cout<<"Mt: "<<Mt.numrows()<<"x"<<Mt.numcols()<<endl;
  946.     // cout<<"y: "<<y.numrows()<<"x"<<y.numcols()<<endl;
  947.     // cout<<"MtM: "<<MtM.numrows()<<"x"<<MtM.numcols()<<endl;
  948.     // cout<<"Mty: "<<MtY.numrows()<<"x"<<MtY.numcols()<<endl;
  949.     return redaug(augment(MtM, MtY),MtM, MtY); //[MtM : MtY]
  950.  
  951. }
  952. void rVals()   //Enter X and Y values for Matrix Regression
  953. {
  954.     cout<<"Enter x values. Type 'd' when finished.\n";
  955.     bool eX = true;
  956.     dataX.resize(row+1, 1);
  957.     int row=0;
  958.     while(eX)
  959.     {
  960.         cout<<">";
  961.  
  962.         string n;
  963.         cin>>n;
  964.         if(n=="d")
  965.         {
  966.             dataX.resize(dataX.numrows()-1,1);
  967.             eX = false;
  968.         }
  969.         else
  970.         {
  971.             double num;
  972.             dataX[row][0]=strtod(n.c_str(),NULL);
  973.             row++;
  974.             dataX.resize(row+1,1);
  975.         }
  976.     }
  977.     //printM(dataX);
  978.     cout<<"Enter y values.\n";
  979.     dataY.resize(row+1, 1);
  980.     row=0;
  981.     for(int i = 0; i<dataX.numrows(); i++)
  982.     {
  983.         cout<<">";
  984.  
  985.         string n;
  986.         cin>>n;
  987.         if(n=="d")
  988.         {
  989.             dataY.resize(dataY.numrows()-1,1);
  990.  
  991.         }
  992.         else
  993.         {
  994.             double num;
  995.             dataY[row][0]=strtod(n.c_str(),NULL);
  996.             row++;
  997.             dataY.resize(row+1,1);
  998.         }
  999.     }
  1000. }
  1001. apmatrix<double> augment(apmatrix<double> A, apmatrix<double> B)   //Instead of inverse for Matrix Regression
  1002. {
  1003.     apmatrix<double>aug;
  1004.  
  1005.     aug = A;
  1006. //A = nxn
  1007. //B = nx1
  1008.     aug.resize(A.numrows(),A.numcols()+B.numcols()); //A augmented with B matrix
  1009. //A.resize(A.numrows(),(2*mtx.numcols())); //Space to add B to A
  1010. //Adds B to A
  1011.     for(int i = 0; i<A.numrows(); i++)
  1012.     {
  1013.         for(int j = A.numcols(); j<aug.numcols(); j++)
  1014.         {
  1015.             aug[i][j]=B[i][j-A.numcols()]; //Adds B to aug
  1016.         }
  1017.  
  1018.     }
  1019.     return aug;
  1020. }
  1021. apmatrix<double> redaug(apmatrix<double> aug, apmatrix<double> A, apmatrix<double> B)   //Reduces Augmented Matrix and gets new inverted matrix on right side
  1022. {
  1023.     apmatrix<double>X;
  1024.     X.resize(B.numrows(), B.numcols());
  1025.     aug = rrE(rE(aug,aug),rE(aug,aug));//Reduces to get X  matrix on right side
  1026.     for(int i = 0; i<aug.numrows(); i++)
  1027.     {
  1028.         for(int j = (A.numcols()); j<aug.numcols(); j++)
  1029.         {
  1030.             X[i][j-(A.numcols())]=aug[i][j]; //Sets X to x values
  1031.         }
  1032.     }
  1033.     return X;
  1034. }
  1035. apmatrix<double> proj(apmatrix<double> v, apmatrix<double> u)   //Projection of u onto v
  1036. {
  1037. //proj(v,u) = [(u dot v)/(v dot v)]*v
  1038.     double a = mult(transpose(u,u),v)[0][0], b = mult(transpose(v,v),v)[0][0], c = (a/b);
  1039.     apmatrix<double>k;
  1040.     k.resize(v.numrows(),v.numcols());
  1041.     for(int i = 0; i<v.numrows(); i++)
  1042.     {
  1043.         for(int j = 0; j<v.numcols(); j++)
  1044.         {
  1045.             k[i][j]=v[i][j]*c;
  1046.         }
  1047.     }
  1048.     return k;
  1049. }
  1050. apmatrix<double> orthogonal(apmatrix<double> A)   //Orthogonal
  1051. {
  1052.     /*
  1053.     n = (number of columns of A)-1
  1054.     x = columns of A
  1055.     v0 = x0
  1056.     v1 = x1 - proj(v0,x1)
  1057.     vn = xn - proj(v0,xn) - proj(v1,xn) - ... - proj(v(n-1),xn)
  1058.     V = {v1, v2, ..., vn} or [v0 v1 ... vn]
  1059.     */
  1060.     apmatrix<double> V, x, v;
  1061.     int n = A.numcols();
  1062.     V.resize(A.numrows(),n);
  1063.     x.resize(A.numrows(), 1);
  1064.     v.resize(A.numrows(),1);
  1065.     for(int i = 0; i<A.numrows(); i++)
  1066.     {
  1067.         x[i][0]=A[i][1];
  1068.         v[i][0]=A[i][0];
  1069.         V[i][0]=A[i][0];
  1070.     }
  1071.     for (int c = 1; c<n; c++)   //Iterates through each col of A as if each was its own matrix
  1072.     {
  1073.         apmatrix<double>vn,vc; //vn = Orthogonalized v (avoiding matrix overwriting of v); vc = previously orthogonalized v
  1074.        // vn.resize(v.numrows(), 1);
  1075.        // for(int l = 0; l<vn.numrows(); l++)
  1076.             vn=x;
  1077.         vc.resize(v.numrows(), 1);
  1078.         for(int i=0; i<c; i++)   //Vn = an-(sigma(t=1, n-1, proj(vt, xn))
  1079.         {
  1080.             for(int k = 0; k<V.numrows(); k++)
  1081.                 vc[k][0] = V[k][i]; //Sets vc to designated v matrix
  1082.             cout<<"c: "<<c<<"\ni: "<<i<<"\nvc x proj\n";
  1083.             apmatrix<double>temp = proj(vc, x);
  1084.             for(int j = 0; j<A.numrows(); j++)
  1085.             {
  1086.                 cout<<vc[j][0]<<" "<<x[j][0]<<" "<<temp[j][0]<<endl;
  1087.                 //vn = add(v,proj(vc,x));
  1088.                 vn[j][0]-=temp[j][0]; //orthogonalize matrix
  1089.             }
  1090.             cout<<endl;
  1091.         }
  1092.         for(int k = 0; k<V.numrows(); k++)
  1093.         {
  1094.             V[k][c]=vn[k][0]; //Subtracts orthogonalized col to V
  1095.             v[k][0]=V[k][c]; //Sets v = new orthogonalized col for
  1096.         }
  1097.         if((c+1)<A.numcols())  //Matrix Out of Bounds Checker
  1098.         {
  1099.             for(int k = 0; k<A.numrows(); k++)
  1100.             {
  1101.                 vn[k][0]=0;
  1102.                 vc[k][0]=0;
  1103.                 x[k][0]=A[k][c+1]; //Moves x onto next v
  1104.             }
  1105.         }
  1106.     }
  1107.     cout<<"\nV:\n";
  1108.     for(int i = 0; i<V.numrows(); i++)
  1109.     {
  1110.         for(int j = 0; j<V.numcols(); j++)
  1111.         {
  1112.             cout<<V[i][j]<<" ";
  1113.         }
  1114.         cout<<endl;
  1115.     }
  1116.     system("PAUSE");
  1117.     return V;
  1118. }
  1119. apmatrix<double> orthonormal(apmatrix<double> V)   //Orthonormal
  1120. {
  1121. // v/||v||
  1122.     apmatrix<double> O = V;
  1123.     for(int c = 0; c<V.numcols(); c++)  //Iterates through each column of V
  1124.     {
  1125.         double mag; //Magnitude of v
  1126.         for(int i = 0; i<O.numrows(); i++)
  1127.         {
  1128.             mag+=pow(O[i][c],2); //sqrt(v1n^2 + v2n^2 + ... + vmn^2)
  1129.         }
  1130.         mag = sqrt(mag);
  1131.         for(int i=0; i<O.numrows(); i++)   //V/||v||
  1132.         {
  1133.             O[i][c]*=(1/mag);
  1134.         }
  1135.     }
  1136.     return O;
  1137. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement