Got LeastSquares working (3ec74f18) · Commits · aflab / astrogeology / ISIS3

isis/src/base/objs/LeastSquares/LeastSquares.cpp

+9 −43

Original line number	Diff line number	Diff line
		@@ -44,6 +44,9 @@ namespace Isis {
		p_sparse = sparse;
		p_sigma0 = 0.;

		p_sparseRows = sparseRows;
		p_sparseCols = sparseCols;


		if (p_sparse) {

		@@ -67,9 +70,6 @@ namespace Isis {
		p_parameterWeights.resize(sparseCols);
		}


		p_sparseRows = sparseRows;
		p_sparseCols = sparseCols;
		}
		p_currentFillRow = -1;
		}
		@@ -156,7 +156,6 @@ namespace Isis {

		p_currentFillRow++;

		// ??? Speed this up using iterator instead of array indices???
		int ncolumns = (int)data.size();

		for(int c = 0; c < ncolumns; c++) {
		@@ -448,13 +447,6 @@ namespace Isis {
		std::cout << p_sparseA.n_rows << ", " << p_sparseA.n_cols << std::endl;
		p_normals = p_sparseA.t()*p_sparseA;

		// // Test for any columns with all 0's
		// // Return column number so caller can determine the appropriate error.
		// arma::umat nonZeros = any(p_normals, 0);
		// for(int c = 0; c < nonZeros.n_rows; c++) {
		// if (nonZeros(0, c) == true) return c + 1;
		// }

		// Create the right-hand-side column vector 'b'
		arma::mat b(p_sparseRows, 1);

		@@ -474,14 +466,15 @@ namespace Isis {
		continue;

		p_normals(i, i) += weight;
		p_ATb(0, i) -= p_epsilonsSparse[i]*weight;
		p_ATb(i, 0) -= p_epsilonsSparse[i]*weight;
		}
		}

		bool status = spsolve(p_xSparse, p_normals, p_ATb, "superlu");

		if (status == false) {
		//TODO what happens if it can't solve?
		QString msg = "Could not solve sparse least squares problem.";
		throw IException(IException::Unknown, msg, _FILEINFO_);
		}

		// Set the coefficients in our basis equation
		@@ -494,23 +487,6 @@ namespace Isis {
		p_epsilonsSparse[i] += p_xSparse[i];
		}

		// test print solution
		// printf("printing solution vector and epsilons\n");
		// for( int a = 0; a < p_sparseCols; a++ )
		// printf("soln[%d]: %lf epsilon[%d]: %lf\n",a,p_xSparse[a],a,p_epsilonsSparse[a]);

		// printf("printing design matrix A\n");
		// for (int i = 0; i < p_sparseRows; i++ )
		// {
		// for (int j = 0; j < p_sparseCols; j++ )
		// {
		// if ( j == p_sparseCols-1 )
		// printf("%20.20e \n",(double)(p_sparseA(i,j)));
		// else
		// printf("%20.20e ",(double)(p_sparseA(i,j)));
		// }
		// }

		// Compute the image coordinate residuals and sum into Sigma0
		// (note this is exactly what was being done before, but with less overhead - I think)
		// ultimately, we should not be using the A matrix but forming the normals
		@@ -546,25 +522,15 @@ namespace Isis {
		p_degreesOfFreedom = p_sparseRows + p_constrainedParameters - p_sparseCols;

		if( p_degreesOfFreedom <= 0.0 ) {
		printf("Observations: %d\nConstrained: %d\nParameters: %d\nDOF: %d\n",
		p_sparseRows,p_constrainedParameters,p_sparseCols,p_degreesOfFreedom);
		p_sigma0 = 1.0;
		}
		else
		else {
		p_sigma0 = p_sigma0/(double)p_degreesOfFreedom;
		}

		// check for p_sigma0 < 0
		p_sigma0 = sqrt(p_sigma0);

		// kle testing - output residuals and some stats
		printf("Sigma0 = %20.10lf\nNumber of Observations = %d\nNumber of Parameters = %d\nNumber of Constrained Parameters = %d\nDOF = %d\n",p_sigma0,p_sparseRows,p_sparseCols,p_constrainedParameters,p_degreesOfFreedom);
		// printf("printing residuals\n");
		// for( int k = 0; k < p_sparseRows; k++ )
		// {
		// printf("%lf %lf\n",p_residuals[k],p_residuals[k+1]);
		// k++;
		// }

		// All done
		p_solved = true;
		return 0;

isis/src/base/objs/LeastSquares/LeastSquares.h

+2 −0

Original line number	Diff line number	Diff line
		@@ -108,6 +108,8 @@ namespace Isis {
		* to reset all solution methods
		* @history 2010-11-22 Debbie A. Cook - Merged with Ken Edmundson version
		* @history 2013-12-29 Jeannie Backer - Improved error messages. Fixes #962.
		* @history 2019-09-05 Makayla Shepherd & Jesse Mapel - Changed sparse solution to use
		* Armadillo library's SuperLU interface instead of GMM.
		*
		*/
		class LeastSquares {

isis/src/base/objs/OverlapNormalization/OverlapNormalization.cpp

+17 −2

Original line number	Diff line number	Diff line
		@@ -177,6 +177,21 @@ namespace Isis {
		throw IException(IException::User, msg, _FILEINFO_);
		}

		if ( method == LeastSquares::SPARSE ) {
		m_offsetLsq = NULL;
		m_gainLsq = NULL;
		delete m_offsetLsq;
		delete m_gainLsq;
		int sparseMatrixRows = m_overlapList.size() + m_idHoldList.size();
		int sparseMatrixCols = m_offsetFunction->Coefficients();
		m_offsetLsq = new LeastSquares(*m_offsetFunction, true, sparseMatrixRows, sparseMatrixCols, true);
		sparseMatrixCols = m_gainFunction->Coefficients();
		m_gainLsq = new LeastSquares(*m_gainFunction, true, sparseMatrixRows, sparseMatrixCols, true);
		const std::vector<double> alphaWeight(sparseMatrixCols, 1/1000.0);
		m_offsetLsq->SetParameterWeights( alphaWeight );
		m_gainLsq->SetParameterWeights( alphaWeight );
		}

		// Calculate offsets
		if (type != Gains && type != GainsWithoutNormalization) {
		// Add knowns to least squares for each overlap
		@@ -252,7 +267,7 @@ namespace Isis {
		m_gainLsq->AddKnown(input, log(tanp), m_weights[overlap]);
		}
		else {
		m_gainLsq->AddKnown(input, 0.0, 1e30); // Set gain to 1.0
		m_gainLsq->AddKnown(input, 0.0, 1e10); // Set gain to 1.0
		}
		}

		@@ -264,7 +279,7 @@ namespace Isis {
		input.resize(m_statsList.size());
		for (int i = 0; i < (int)input.size(); i++) input[i] = 0.0;
		input[hold] = 1.0;
		m_gainLsq->AddKnown(input, 0.0, 1e30);
		m_gainLsq->AddKnown(input, 0.0, 1e10);
		}

		// Solve the least squares and get the gain coefficients to apply to the

isis/src/base/objs/OverlapNormalization/OverlapNormalization.h

+2 −0

Original line number	Diff line number	Diff line
		@@ -70,6 +70,8 @@ namespace Isis {
		* message. Fixes #962,
		* @history 2013-02-14 Steven Lambright - Added SolutionType GainsWithoutNormalization. Fixes
		* #911.
		* @history 2019-09-05 Makayla Shepherd & Jesse Mapel - Changed weight for hold images from
		* 1E30 to 1E10 to avoid poorly conditioned normal matrix.
		*/

		class OverlapNormalization {