23 #include "jama/jama_svd.h"
24 #include "jama/jama_qr.h"
27 #include "gmm/gmm_superlu_interface.h"
42 int sparseRows,
int sparseCols,
bool jigsaw) {
56 if (sparseRows == 0 || sparseCols == 0) {
57 QString msg =
"If solving using sparse matrices, you must enter the "
58 "number of rows/columns";
63 gmm::resize(
p_sparseA, sparseRows, sparseCols);
64 gmm::resize(
p_normals, sparseCols, sparseCols);
65 gmm::resize(
p_ATb, sparseCols, 1);
76 p_sparseRows = sparseRows;
77 p_sparseCols = sparseCols;
79 p_currentFillRow = -1;
120 QString msg =
"Number of elements in data does not match basis [" +
168 int ncolumns = (int)data.size();
170 for(
int c = 0; c < ncolumns; c++) {
185 if((row >=
Rows()) || (row < 0)) {
186 QString msg =
"Index out of bounds [Given = " +
toString(row) +
"]";
200 if((row >=
Rows()) || (row < 0)) {
201 QString msg =
"Index out of bounds [Given = " +
toString(row) +
"]";
237 if((method ==
SPARSE && p_sparseRows == 0) ||
240 QString msg =
"No solution available because no input data was provided";
247 else if(method ==
QRD) {
250 else if(method ==
SPARSE) {
272 for(
int r = 0; r < A.dim1(); r++) {
274 for(
int c = 0; c < A.dim2(); c++) {
285 JAMA::SVD<double> svd(A);
287 TNT::Array2D<double> V;
291 TNT::Array2D<double> invS;
294 for(
int i = 0; i < invS.dim1(); i++) {
295 if(invS[i][i] != 0.0) invS[i][i] = 1.0 / invS[i][i];
299 TNT::Array2D<double> U;
301 TNT::Array2D<double> transU(U.dim2(), U.dim1());
303 for(
int r = 0; r < U.dim1(); r++) {
304 for(
int c = 0; c < U.dim2(); c++) {
305 transU[c][r] = U[r][c];
310 TNT::Array2D<double> VinvS = TNT::matmult(V, invS);
311 TNT::Array2D<double> Aplus = TNT::matmult(VinvS, transU);
316 for(
int r = 0; r < (int)
p_expected.size(); r++) {
320 TNT::Array2D<double> coefs = TNT::matmult(Aplus, b);
325 QString msg =
"Unable to solve least-squares using SVD method. No "
326 "solution available. Not enough knowns or knowns are "
327 "co-linear ... [Unknowns = "
334 std::vector<double> bcoefs;
335 for (
int i = 0; i < coefs.dim1(); i++) bcoefs.push_back(coefs[i][0]);
340 for(
int i = 0; i < (int)
p_input.size(); i++) {
375 for(
int r = 0; r < A.dim1(); r++) {
377 for(
int c = 0; c < A.dim2(); c++) {
387 JAMA::QR<double> qr(A);
391 for(
int r = 0; r < (int)
p_expected.size(); r++) {
398 int full = qr.isFullRank();
400 QString msg =
"Unable to solve-least squares using QR Decomposition. "
401 "The upper triangular R matrix is not full rank";
405 TNT::Array1D<double> coefs = qr.solve(b);
408 std::vector<double> bcoefs;
409 for(
int i = 0; i < coefs.dim1(); i++) {
410 bcoefs.push_back(coefs[i]);
415 for(
int i = 0; i < (int)
p_input.size(); i++) {
469 for(
int c = 0; c < p_sparseCols; c++) {
470 numNonZeros = gmm::nnz(gmm::sub_matrix(
p_normals,
471 gmm::sub_interval(0,p_sparseCols),
472 gmm::sub_interval(c,1)));
474 if(numNonZeros == 0)
return c + 1;
478 gmm::dense_matrix<double> b(p_sparseRows, 1);
481 for (
int r = 0; r < p_sparseRows; r++ )
489 for(
int i = 0; i < p_sparseCols; i++) {
518 for(
int i = 0; i < p_sparseCols; i++ )
548 for (
int i = 0; i < p_sparseRows; i++ ) {
557 double constrained_vTPv = 0.0;
559 for (
int i = 0; i < p_sparseCols; i++ ) {
574 printf(
"Observations: %d\nConstrained: %d\nParameters: %d\nDOF: %d\n",
658 for (
int i = 0; i < p_sparseCols; i++ )
701 void LeastSquares::Reset ()
707 p_currentFillRow = -1;
734 QString msg =
"Unable to evaluate until a solution has been computed";
751 QString msg =
"Unable to return residuals until a solution has been computed";
771 QString msg =
"Unable to return residuals until a solution has been computed";