import util::Math;
import Prelude;
import vis::Figure;
import vis::Render;

public rel[real,real,real] QRdecomposition(rel[real x, real y, real v] matrix){
	//orthogonalcolumns
	oc = domainR(matrix, {0.0});
	for (x <- sort(toList(domain(matrix)-{0.0}))){
		c = domainR(matrix, {x});
		o = domainR(oc, {x-1});
		
		for (n <- [1.0 .. x]){
			o = domainR(oc, {n-1});
			c = matrixSubtract(c, matrixMultiplybyN(o, matrixDotproduct(o, c)/matrixDotproduct(o, o)));
			}
			
		oc += c;
	}
	
	Q = {};
	//from orthogonal to orthonormal columns
	for (el <- oc){
		c = domainR(oc, {el[0]});
		Q += matrixNormalize({el}, c);
	}
	
	//from Q to R
	R= matrixMultiplication(matrixTranspose(Q), matrix);
	R= {<x,y,toReal(round(v))> | <x,y,v> <- R};
	
	println("Q:");
	iprintlnExp(Q);
	println();
	println("R:");
	return R;
}

//a function that takes the transpose of a matrix, see also Rosetta Code problem "Matrix transposition"
public rel[real, real, real] matrixTranspose(rel[real x, real y, real v] matrix){
	return {<y, x, v> | <x, y, v> <- matrix};
}

//a function to normalize an element of a matrix by the normalization of a column
public rel[real,real,real] matrixNormalize(rel[real x, real y, real v] element, rel[real x, real y, real v] column){
	normalized = 1.0/nroot((0.0 | it + v*v | <x,y,v> <- column), 2);
	return matrixMultiplybyN(element, normalized);
}

//a function that takes the dot product, see also Rosetta Code problem "Dot product"
public real matrixDotproduct(rel[real x, real y, real v] column1, rel[real x, real y, real v] column2){
	return (0.0 | it + v1*v2 | <x1,y1,v1> <- column1, <x2,y2,v2> <- column2, y1==y2);
}

//a function to subtract two columns
public rel[real,real,real] matrixSubtract(rel[real x, real y, real v] column1, rel[real x, real y, real v] column2){
	return {<x1,y1,v1-v2> | <x1,y1,v1> <- column1, <x2,y2,v2> <- column2, y1==y2};
}

//a function to multiply a column by a number
public rel[real,real,real] matrixMultiplybyN(rel[real x, real y, real v] column, real n){
	return {<x,y,v*n> | <x,y,v> <- column};
}

//a function to perform matrix multiplication, see also Rosetta Code problem "Matrix multiplication".
public rel[real, real, real] matrixMultiplication(rel[real x, real y, real v] matrix1, rel[real x, real y, real v] matrix2){
	if (max(matrix1.x) == max(matrix2.y)){
		p = {<x1,y1,x2,y2, v1*v2> | <x1,y1,v1> <- matrix1, <x2,y2,v2> <- matrix2};

		result = {};
		for (y <- matrix1.y){
			for (x <- matrix2.x){
				v = (0.0 | it + v | <x1, y1, x2, y2, v> <- p,  x==x2 && y==y1, x1==y2 && y2==x1);
				result += <x,y,v>;
			}
		}
		return result;
	}
	else throw "Matrix sizes do not match.";
} 	

// a function to visualize the result
public void displayMatrix(rel[real x, real y, real v] matrix){
	points = [box(text("<v>"), align(0.3333*(x+1),0.3333*(y+1)),shrink(0.25)) | <x,y,v> <- matrix];
	render(overlay([*points], aspectRatio(1.0)));
}

//a matrix, given by a relation of <x-coordinate, y-coordinate, value>.
public rel[real x, real y, real v] matrixA = {
<0.0,0.0,12.0>, <0.0,1.0, 6.0>, <0.0,2.0,-4.0>,
<1.0,0.0,-51.0>, <1.0,1.0,167.0>, <1.0,2.0,24.0>,
<2.0,0.0,4.0>, <2.0,1.0,-68.0>, <2.0,2.0,-41.0>
};