From 27d7f9dac270da1de1a6d17dee05b35f1c7c7560 Mon Sep 17 00:00:00 2001
From: nrubin29 <nrubin29@gmail.com>
Date: Sat, 23 Dec 2017 22:38:56 -0500
Subject: [PATCH] Added QR factorization.

---
 common.py   |  1 +
 rules.py    |  6 +++---
 vartypes.py | 61 ++++++++++++++++++++++++++++++++++++++++++++++++++---
 3 files changed, 62 insertions(+), 6 deletions(-)

diff --git a/common.py b/common.py
index 09b39fb..42ee269 100644
--- a/common.py
+++ b/common.py
@@ -48,6 +48,7 @@ def _str(self, node, depth=0) -> str:
     (r'solve',              'OPR'),
     (r'eval',               'OPR'),
     (r'ls',                 'OPR'),
+    (r'qr',                 'OPR'),
     (r'[a-zA-Z_]+',         'IDT'),
     (r'=',                  'EQL'),
     (r'\+',                 'ADD'),
diff --git a/rules.py b/rules.py
index 704af4d..321a043 100644
--- a/rules.py
+++ b/rules.py
@@ -4,7 +4,7 @@
 from typing import List
 
 from common import Token
-from vartypes import VariableValue, NumberValue, MatrixRowValue, MatrixValue, Value, OperatorBodyValue
+from vartypes import VariableValue, NumberValue, MatrixRowValue, MatrixValue, Value, TupleValue
 
 
 def flatten(l):
@@ -35,8 +35,8 @@ def mbd(values: List[Value], _) -> MatrixValue:
     return MatrixValue(values)
 
 
-def opb(values: List[Value], _) -> OperatorBodyValue:
-    return OperatorBodyValue(values)
+def opb(values: List[Value], _) -> TupleValue:
+    return TupleValue(values)
 
 
 def add(operands: List[Value], operator: Token) -> Value:
diff --git a/vartypes.py b/vartypes.py
index dd8ad7a..2248ef3 100644
--- a/vartypes.py
+++ b/vartypes.py
@@ -107,6 +107,9 @@ def exp(self):
     def identity(self):
         return MatrixValue([[1 if col is row else 0 for col in range(int(self.value))] for row in range(int(self.value))])
 
+    def zeroes(self, other):
+        return MatrixValue([[0 for _ in range(int(other.value))] for _ in range(int(self.value))])
+
 
 class MatrixValue(Value):
     def __init__(self, data):
@@ -123,13 +126,25 @@ def __init__(self, data):
     def __str__(self):
         return '[\n' + '\n'.join(['[' + ', '.join(map(lambda cell: str(round(cell, 5)), row)) + ']' for row in self.value]) + '\n]'
 
+    def sub(self, other):
+        if isinstance(other, MatrixValue):
+            if len(self.value) != len(other.value) or len(self.value[0]) != len(other.value[0]):
+                raise EvaluationException('Attempted to subtract two matrices of different dimensions')
+
+            return MatrixValue([[self.value[row][col] - other.value[row][col] for col in range(len(self.value[row]))] for row in range(len(self.value))])
+
+        raise EvaluationException('Cannot sub {} and {}'.format(self.type, other.type))
+
     def mul(self, other):
         if isinstance(other, NumberValue):
-            # Number * Matrix
+            # Matrix * Number
             return MatrixValue([[cell * other.value for cell in row] for row in self.value])
 
         elif isinstance(other, MatrixValue):
             # Matrix * Matrix
+            if len(self.value[0]) != len(other.value):
+                raise EvaluationException('Cannot multiply matrices of dimensions {} and {}'.format((len(self.value), len(self.value[0])), (len(other.value), len(other.value[0]))))
+
             result = [[0 for _ in range(len(other.value[0]))] for _ in range(len(self.value))]
 
             for i in range(len(self.value)):
@@ -141,7 +156,15 @@ def mul(self, other):
 
         else:
             raise EvaluationException('Cannot mul {} and {}'.format(self.type, other.type))
-    
+
+    def div(self, other):
+        if isinstance(other, NumberValue):
+            # Matrix / Number
+            return MatrixValue([[cell / other.value for cell in row] for row in self.value])
+
+        else:
+            raise EvaluationException('Cannot div {} and {}'.format(self.type, other.type))
+
     def det(self):
         return NumberValue(self._det(self.value))
 
@@ -207,6 +230,38 @@ def solve(self, other):
     def ls(self, other):
         return (self.trans().mul(self)).inv().mul(self.trans()).mul(other)
 
+    def norm(self):
+        return NumberValue(math.sqrt(sum([sum([col * col for col in row]) for row in self.value])))
+
+    def squeeze(self):
+        return MatrixValue([[cell for row in self.value for cell in row]])
+
+    def _col(self, col):
+        """ Isolates an individual column from the matrix """
+        return MatrixValue([[row[col]] for row in self.value])
+
+    def qr(self):
+        m = NumberValue(len(self.value))
+        n = NumberValue(len(self.value[0]))
+
+        Q = m.zeroes(m).value
+        R = n.zeroes(n).value
+
+        for j in range(n.value):
+            v = self._col(j)
+
+            for i in range(j):
+                R[i][j] = MatrixValue(Q)._col(i).trans().mul(self._col(j)).value[0][0]
+                v = v.squeeze().sub(MatrixValue(Q)._col(i).trans().mul(NumberValue(R[i][j])))
+
+            R[j][j] = v.norm().value
+
+            val = v.div(NumberValue(R[j][j])).squeeze().value[0]
+            for row in range(len(Q)):
+                Q[row][j] = val[row]
+
+        return TupleValue([Q, R])
+
 
 class MatrixRowValue(Value):
     def __init__(self, data):
@@ -219,7 +274,7 @@ def __init__(self, data):
             self.value = data
 
 
-class OperatorBodyValue(Value):
+class TupleValue(Value):
     def __init__(self, args: List):
         super().__init__()