94 lines
3.1 KiB
Plaintext
94 lines
3.1 KiB
Plaintext
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>
|
|
};
|